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SUMMARY 

Both  mortality  rate  and  radial  growth  of  high 
elevation  (>900  m)  red  spruce-Fraser  fir  forests 
of  the  southern  Appalachians  have  experienced 
change  since  approximately  1960.  Scientific 
interest  in  a  study  of  these  forests  have 
increased  because  atmospheric  pollution  is  a 
possible  cause  of  the  change.  Scientists  with 
statistical  and  biological  expertise 
independently  analyzed  a  tree  ring  data  set 
collected  by  the  Tennessee  Valley  Authority  and 
the  National  Park  Service.  The  objective  of  the 
analysis  was  to  develop  new  or  improved 
techniques  for  extracting  information  from  such 
data;  tree  rings  represent  a  natural  data  storage 
system  that  is  one  of  the  few  sources  of  long- 
term  information  for  these  forests. 

Although  no  definite  statements  are  made  about 
the  role  of  atmospheric  deposition  in  observed 
forest  decline,  the  results  should  contribute  to 
the  success  of  future  research.  The  four 
techniques  employed  in  the  study  involved:  (1)  a 
dendrochronological  approach  employing  spline 
detrending  and  multiple  regression  to  study  the 
effects  of  climate  on  ring  width,  (2)  an 
application  of  fractals  to  study  the  dependence 
of  variance  on  mean  ring  width  over  time,  (3)  an 
approach  that  combined  Box-Jenkins  methods  and 
spatial  analysis,  and  (4)  a  method  of  studying 
time  dependence  of  ring  width  on  climate  using 
the  Kalman  filter. 
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INTRODUCTION 


The  Institute  for  Quantitative  Studies  at  the 
USDA  Forest  Service,  Southern  Forest  Experiment 
Station,  has  been  engaged  in  a  study  of 
statistical  methods  used  in  researching 
Atmospheric  Deposition  Influences  on  Forests 
(ADIF) .  The  study  was  funded  by  the  National 
Vegetation  Survey,  which  is  under  the  Natio..  1 
Acid  Precipitation  Assessment  Program  (NAPAP) . 
The  study  began  with  a  distributed  seminar;  for 
10  weeks  articles  about  ADIF  were  sent  to  a 
number  of  participants  who  returned  comments  each 
week.  A  final  report  was  produced  (Kiester  and 
others  1985)  consisting  of  critiques  of  past  ADIF 
studies,  suggestions  for  additional  reading,  and 
philosophy  about  the  type  and  quality  of  research 
needed  in  the  ADIF  area. 

The  study  indicated  that  tree  ring  analyses 
held  promise  for  the  study  of  ADIF,  but  further 
development  of  appropriate  statistical  methods 
would  be  useful.  Therefore  the  idea  of 
"replicated-statisticians"  was  employed.  Because 
there  are  several  ways  to  approach  a  tree  ring 
analysis,  it  was  probable  that  allowing  a  number 
of  individuals  to  work  independently  would  yield 
some  interesting  new  dendrochronological 
techniques.  Red  spruce  (Picea  rubens  Sarg.)  was 
chosen  for  this  study  because  claims  were 
previously  made  that  spruce  forests  were 
experiencing  unexplained  growth  declines  in  both 
the  Northeast  and  the  South  (Hornbeck  and  Smith 
1985;  Adams  and  others  1985).  Although  the  study 
was  intended  to  develop  statistical  methods  and 
not  explanations  of  growth  declines,  using  this 
data  allowed  for  that  possibility 

This  study  was  a  joint  endeavor  by  the  USDA 
Forest  Service,  the  National  Park  Service  (NPS), 
and  the  Tennessee  Valley  Authority  (TVA) . 
Funding  came  from  the  Forest  Service  and  the  TVA, 
data  from  the  TVA  and  NPS,  and  statistical 
analyses  were  conducted  by  the  Forest  Service. 
Agreements  were  made  with  the  following 
scientists  to  perform  analyses  and  provide 
individual  reports : 

(1)  Edward  R.  Cook,  Ph.D. 
Tree -Ring  Laboratory 

Lamont-Doherty  Geological  Observatory 
of  Columbia  University 

Palisades,  New  York  10964 

(2)  Keith  Ord,  Ph.D.  and  Janice  Derr,  Ph.D. 
Department  of  Management  Science 
Pennsylvania  State  University 
University  Park,  PA  16802 

(3)  Robin  A.  J.  Taylor,  Ph.D. 
Department  of  Entomology 
106  Patterson  Building 
Pennsylvania  State  University 
University  Park,  PA  16802 

(4)  Paul  C.  Van  Deusen,  Ph.D. 
Institute  for  Quantitative  Studies 
Southern  Forest  Experiment  Station 
701  Loyola  Avenue 

New  Orleans,  LA  70113 

Dr.  Cook  is  considered  to  be  one  of  the  leaders 
in  the  field  of  dendrochronology,  and  his 
analysis  therefore  includes  the  most  currently 
accepted  techniques.   He  concluded  that  there  is 


evidence  of  anomalous  behavior  in  the  red  spruce 
forests  of  the  Great  Smoky  Mountains  and  suggests 
that  this  is  partly  due  to  warmer  summer 
temperatures  in  recent  years.  He  also  states 
that  the  situation  in  southern  red  spruce  is  not 
similiar  to  that  in  northern  red  spruce. 

Drs .  Ord  and  Derr  performed  a  spatial  analysis 
on  the  data.  They  concluded  that  there  was  a 
tendency  for  ring  widths  to  diminish  in  recent 
years  and  that  there  is  a  strong  spatial 
dependence  in  forecast  residuals  that  cannot  be 
explained  by  geographic  or  biotic  factors  alone. 
They  did  not  consider  climate,  which  they  mention 
could  explain  some  of  the  remaining  spatial 
dependence.  The  analysis  performed  here  is  novel 
for  the  field  of  dendrochronology  and  may  lead  to 
useful  results  in  the  future. 

To  study  the  dependence  of  variance  on  the  mean 
of  the  data,  Dr.  Taylor  investigated  the  use  of 
"fractals,"  a  term  coined  to  denote  fractional 
dependence.  He  concluded  that  the  change  in 
fractional  dimension  over  time  may  be  due  to 
successional ,  climatic,  or  anthropogenic 
influences.  This  technique  has  not  been  applied 
previously  to  tree  ring  data  and  may  show  promise 
after  further  development. 

Dr.  Van  Deusen  analyzed  the  data  using  the 
Kalman  filter  technique,  which  is  commonly  used 
in  engineering  applications.  At  the  time  of  this 
study,  the  method  had  not  been  used  in 
dendrochronology,  although  scientists  in  the 
Netherlands  (Visser  1986)  have  recently  published 
a  paper  on  independent  applications  of  the  method 

to  tree  rings.  The  Kalman  Filter  allowed  the 
climatic  data  to  be  modeled  dynamically  so  that 
its  effect  over  time  could  be  studied.  He 
concluded  that  these  trees  have  become 
progressively  more  sensitive  to  climate  since  the 
late  1950' s.  This  increased  sensitivity  may 
coincide  with  insect-caused  thinning  in  the 
stands . 

All  the  studies  concluded  that  the  growth 
patterns  in  these  stands  have  changed  in  the  last 
20  years.  The  causes  of  these  changes  are 
uncertain,  but  the  sensitivity  of  the  stand  to 
climate  appears  to  be  increasing.  To  enhance  the 
reader's  ability  to  interpret  the  various 
analyses ,  the  individual  reports  are  prefaced  by 
a  review  by  Elizabeth  Groton  and  Christopher 
Eagar  of  the  geographical  and  biological 
background  of  the  Southern  Appalachian  Spruce  Fir 
Forest. 

References  cited  in  this  section  are:  Adams, 
H.S.;  Stephenson,  S.L.;  Biasing,  T.J.;  Duvick, 
D.N.  1985.  Growth- trend  declines  of  spruce  and 
fir  in  mid-Appalachian  subalpine  forest. 
Environmental  and  Experimental  Botany.  25(4): 
315-325.  Hornbeck,  J.W.;  Smith,  R.B.  1985. 
Documentation  of  red  spruce  growth  decline. 
Canadian  Journal  of  Forest  Research.  15:  1199- 
1201.  Kiester,  A.R. ;  Van  Deusen,  P.C.;  Dell, 
T.R.  1985.  Status  of  the  concepts  and 
methodlogies  used  in  the  study  of  the  effects  of 
atmospheric  deposition  on  forests.  Internal 
report  for  the  National  Vegetation  Survey  of 
NAPAP.  Visser,  H.  1986.  Analysis  of  tree  ring 
data  using  the  Kalman  Filter  techniques.  IAWA 
Bulletin  n.s.,  Vol.  7(4)  289-297.  Published  at 
the  Rijksherbarium,  Netherlands. 


Southern  Appalachian  Red  Spruce- -Fraser  Fir  Forests 

Elizabeth  Groton  and  Christopher  Eagar 


GEOGRAPHICAL  AND  BIOLOGICAL 
BACKGROUND 

The  southern  Appalachian  red  spruce -Fraser  fir+ 
forests  are  found  in  southwestern  Virginia, 
western  North  Carolina,  and  eastern  Tennessee 
(fig.  1).  Total  area  of  these  southern 
Appalachian  spruce -fir  forests  is  estimated  at 
26,577  hectares,  with  19,755  hectares  occurring 
within  the  Great  Smoky  Mountains  National  Park 
(GSMNP).  Extensive  logging  and  other 
disturbances  early  in  the  20th  century  have 
reduced  the  extent  of  the  southern  Appalachian 
spruce-fir  forests.  Today  this  forest  type 
occurs  on  high-elevation  peaks  (above  1,370  m)  in 
islandlike  patches. 

Species'   composition   in   the   southern 
Appalachian  forests  changes  with  the  elevational 
gradient.   At  lower  elevations  (1,370-1,580  m)  in 
undisturbed  forests  such  as  those  of  the  Great 
Smoky  Mountains,  red  spruce  (Picea  rub ens  Sarg.) 
is  found  in  combination  with  northern  hardwood 
type  species:   maple  (Acer  spp . ) ,  American  beech 
(Fa pus  grandifolia  Ehrh.),  yellow  birch  (Betula 
alleghaniensis  Britton) ,  eastern  hemlock  (Tsuga 
candensis  (L.)  Carr.,  northern  red  oak  (Ouercus 
rubra  L.),  Carolina  silverbell  (Halesia  Carolina 
L.),   and   yellow   buckeye   (Aesculus   octandra 
Marsh.).     As   elevation   increases,   the   fir 
component  gains  importance,  and  forest 
composition  changes  to  predominantly  red  spruce  - 
Fraser  fir.   Red  spruce  occurs  less  frequently  at 
the  highest  end  of  the  gradient  (above  1,890  m) , 
giving  way  to  essentially  pure  Fraser  fir  (Abies 
fraseri  (Pursh)  Poir.)  stands  on  mountain  tops 
(Whittaker  1956) .    Mountains  that  were  logged 
during, the  early  part  of  this  century  and  did  not 
experience  postlogging  slash  fires  are  dominated 
by  Fraser  fir.   This  includes  most  of  the  Black 
Mountains,  Balsam  Mountains,  Roan  Mountain,  and 
Mount  Rogers . 

Red  spruce  grows  larger  and  lives  longer  than 
Fraser  fir,  but  Fraser  fir  grows  more  rapidly  and 
produces  more  prolific  seed  crops  than  red 
spruce.  Red  spruce  can  live  for  over  350  years, 
grow  to  40  meters  in  height,  and  have  diameters 
at  breast  height  (d.b.h.)  in  excess  of  1  meter. 
Fraser  fir  seldom  lives  longer  than  150  years  and 
attains  a  maximum  height  of  25  meters  and  d.b.h. 
of  50  centimeters.  Oosting  and  Billings  (1951) 
found  five  times  more  Fraser  fir  than  red  spruce 
seedlings  in  old- growth  stands  in  the  Great  Smoky 
Mountains .  Both  species  are  extremely  shade 
tolerant  and  are  capable  of  resuming  normal 
growth  after  50  years  of  suppression. 


"'Tied  spruce-Frase.r  fir  forests  will  generally 
be  referred  to  as  simply  spruce-fir  forests. 


THE  BALSAM  WOOLY  ADELGID 

Most  spruce- fir  forests  of  the  southern 
Appalachians  have  been  recently  disturbed  by  the 
extensive  mortality  of  Fraser  fir  caused  by  an 
introduced  insect,  the  balsam  woolly  adelgid 
(Adelges  piceae) .  This  pest,  a  native  of  Europe, 
is  a  tiny  insect  that  feeds  on  the  bark  of  true 
firs  (Abies  spp.).  Fraser  fir  is  quickly  killed 
by  the  balsam  woolly  adelgid.  Mortality  occurs 
between  3  and  9  years  from  the  time  of  initial 
infestation,  depending  on  the  size  and  vigor  of 
the  tree  (Amman  and  Speers  1965) .  Tree  death  is 
caused  by  the  diffusion  of  compounds  secreted  by 
the  adelgid  into  the  xylem  during  feeding,  which 
causes  formation  of  premature  heartwood. 
Translocation  of  water  and  minerals  to  the  crown 
are  greatly  reduced,  causing  water  stress  and 
eventual  death  of  the  tree  (Puritch  1971,  1973, 
1977,  Puritch  and  Johnson  1971,  Puritch  and  Petty 
1971). 

The  balsam  woolly  adelgid  was  first  identified 
in  North  America  in  1908  on  balsam  fir  (Abies 
balsamea)  in  Maine  (Kotinsky  1916) .  The  adelgid 
has  caused  extensive  mortality  to  balsam  fir 
throughout  eastern  Canada;  however,  infestations 
have  not  progressed  more  than  80  kilometers 
inland  from  the  coast  because  of  extreme  inland 
winter  conditions  (Balch  1952,  Schooley  and 
Bryant  1978) .  The  balsam  woolly  adelgid  is  not 
present  in  the  northern  Appalachian  spruce -fir 
forest . 

The  balsam  woolly  adelgid  was  detected  in  the 
southern  Appalachians  on  Mount  Mitchell,  North 
Carolina,  in  1957  (Speers  1958).  Subsequent 
surveys  revealed  that  the  adelgid  had  spread 
throughout  the  entire  3,035  hectares  of  Fraser 
fir  type  in  the  Black  Mountains  (Nagel  1959). 
High  mortality  of  Fraser  fir  and  widespread 
adelgid  distribution  indicated  establishment 
prior  to  1957,  perhaps  as  early  as  1940.  Balsam 
woolly  adelgids  were  detected  in  1962  on  Roan 
Mountain  (Ciesla  and  Buchanan  1962) ,  and  in  1963 
infestations  were  located  on  Grandfather  Mountain 
and  on  Mount  Sterling  in  the  GSMNP. 

The  adelgid  arrived  in  the  Clingman's  Dome  area 
of  the  GSMNP  in  the  early  1970' s,  and  tree 
mortality  began  there  in  the  late  1970' s. 
Surveys  found  the  adelgid  in  the  Balsam  Mountains 
and  the  nearby  Plott  Balsams  of  North  Carolina  in 
1968  (Rauschenberger  and  Lambert  1970) .  The 
balsam  woolly  adelgid  was  not  found  on  Mount 
Rogers,  Virginia,  until  1979;  however,  subsequent 
stem  analysis  of  several  trees  within  the 
infested  areas  revealed  adelgid-caused  red  wood 
beginning  in  1962  (Lambert  and  others  1980)  in 
the  annual  rings. 

By  1984  and  1985,  the  balsam  woolly  adelgid  had 
caused  extensive  damage  throughout  the  Black 
Mountains.   Balsam   Mountains,   Plott   Balsams, 


Elizabeth  Groton  is  forest  biometrician  for  the  Tennessee  Valley  Authority, 
Norris,  TN.  Christopher  Eagar  is  ecologist  for  the  National  Park  Service, 
Gatlinburg,  TN. 
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Figure  I. --Area    covered   by   the   southern  Appalachian   red   spruce- -Fraser  fir  forests 


Grandfather  Mountain,  and  most  of  the  Great  Smoky 
Mountains.  Limited  use  of  insecticides  at  Roan 
Mountain  reduced  fir  mortality  in  accessible 
areas,  although  nontreated  areas  experienced 
heavy  damage.  In  the  Clingman's  Dome  area  of  the 
GSMNP,  adelgid  infestations  had  caused 
significant  fir  mortality  at  elevations  below 
1,830  meters  and  minimal  damage  above  this 
elevation.  Mount  Rogers  had  suffered  the  least 
Fraser  fir  mortality  of  the  southern  Appalachian 
spruce-fir  forests.  There  were  isolated,  dead 
Fraser  fir  in  areas  known  to  have  been  infested 
for  23  years,  but  even  within  these  areas  the 
impact  of  the  adelgid  was  surprisingly  low. 
Possible  explanations  for  this  anomalous 
condition  on  Mount  Rogers  include:  a  genetic 
based  difference  in  defense  mechanism  to  adelgid 
infestation  of  this  fir  population,  a  reduction 
in  the  toxicity  of  the  secretions  of  the  Mount 
Rogers  adelgid  population,  or  a  combination  of 
both  possibilities. 


RESEARCH  AND  STUDY  ENDEAVORS 

National  Park  Service  and 
Tennessee  Valley  Authority 

Increased  mortality  (Siccama  and  others  1982, 
Scott  and  others  1984,  Vogelmann  and  others  1985) 
and  apparent  reductions  in  radial  increment 
(Adams  and  others  1985,  Bruck  1986,  McLaughlin 
and  others  1983)  in  the  high-elevation  spruce-fir 
forests  of  the  Eastern  United  States  prompted  two 
studies  in  the  southern  Appalachians.  The 
studies  were  designed  to  assess  the  current 
condition  of  these  forests,  relate  observed 
decline  symptoms  to  site  characteristics,  and 
provide  baseline  data  to  monitor  future  changes 
in  the  forest  condition. 

The  first  study,  conducted  by  the  National  Park 
Service  (NPS) ,  began  in  the  summer  of  1984  in  the 


spruce -fir  forests  of  the  GSMNP.  The  second 
study,  conducted  by  the  Tennessee  Valley 
Authority  (TVA)  in  the  autumn  of  1984, 
established  plots  throughout  the  range  of  the 
southern  Appalachian  spruce -fir  type,  excluding 
the  GSMNP.  Intending  to  combine  data  sets  for 
future  analyses,  both  agencies  collaborated  on 
sampling  design  in  order  to  ensure  that  similar 
data  were  collected. 

Problems  of  assessing  change  in  forest 
productivity  prompted  the  TVA  and  NPS  to 
establish  permanent  vegetation  plots  in  the 
southern  range  of  the  spruce -fir  type.  Plot 
establishment  was  also  influenced  by  the  need  for 
additional  information  on  stand  dynamics  in  the 
spruce-fir  forests.  This  led  to  a  collaborative 
study  between  the  TVA,  NPS,  and  the  Forest 
Service.  Data  collected  by  the  TVA  and  NPS 
included  tree  increment  core  data  and  detailed 
plot  information.  The  tree  core  data  and  plot 
data  were  summarized  and  made  available  to  a 
number  of  individual  researchers  for  independent 
analysis.  Plot  information  included  elevation, 
latitude  and  longitude,  live  and  dead  basal  area, 
and  stand  density.  Regional  climatic  data 
(monthly  averages  of  precipitation  and 
temperature  since  1933)  were  also  included  in  the 
data. 

Sampling  procedures  utilized  by  the  NPS  and  TVA 
were  basically  the  same  (table  1).  Both  agencies 
used  stratified  random  sampling,  locating  plots 
on  aerial  photographs  and  topographic  maps .  Data 
that  were  collected  from  the  plots  included  site 
characteristic  data  such  as  slope,  aspect, 
topographic  location,  and  descriptions  of 
understory  vegetation.  Individual  trees  were 
mapped,  measured,  and  evaluated  for  decline 
symptomology .  Quantitative  assessments  were 
made  of  mortality  and  regeneration.  Increment 
cores  were  collected  from  five  dominant  or 
codominant  trees  at  each  site.  Two  cores  per 
tree  were  taken  at  d.b.h. 


Table  1 Comparison    of   Tennessee    Valley  Authority    (TVA)    and   National 

Park  Service    (NPS)    spruce- fir  sampling  procedures 


Variable 


TVA 


NPS 


Plot  Location 


Stratified  random 
(Strata:  Elevation 
and  dry/wet) 


Stratified  random 
(Strata:  Elevation, 
topo  position, 
macro-aspect) 


Elevation 


Variable  within 
strata 


Held  constant  at 
a  set  strata 


Plot  Size 


Circular,  0.08  ha 


Square,  0.04  ha 


Overstory  stand 
data+ 

Site  charater- 
istic  data+ 


Increment  cores 


Cores  taken  from  five 
dominant  or  codominant 
trees  from  outside  the 
plot  and  extended  to 
tree  center. 


Cores  taken  from  five 
dominant  trees  from 
within  plot,  cores 
not  extended  to  tree 
center . 


Overstory  stand  data  and  site  characteristic  data  were  basically 


Identical  for  both  TVA  and  NPS. 


The  results  of  the  analyses  of  tree  core  data 
may  provide  insight  into  the  question  of  whether 
or  not  the  southern  Appalachian  red  spruce  and 
Fraser  fir  are  experiencing  a  decline  that  cannot 
be  attributed  to  natural  stresses.  The  plots 
established  by  the  NPS  and  TVA  will  continue  to 
be  remeasured  in  order  to  monitor  future  changes 
in  stand  productivity. 

Other  Studies 

Several  studies  have  used  annual  radial 
increments  from  tree  cores  to  evaluate  changes  in 
the  growth  rate  of  red  spruce  and  Fraser  fir  from 
several  sites  in  the  southern  Appalachians. 
These  studies  indicated  an  abrupt  shift  to  narrow 
growth  rings  beginning  in  the  late  1960 's  to 
early  1970 's  (Adams  and  others  1985,  Bruck  1986, 
McLaughlin  and  others  1983)  for  red  spruce  and, 
to  a  lesser  extent,  for  Fraser  fir.  This 
tendency  was  more  drastic  at  high  elevations 
(Adams  and  others  1985) .  The  annual  growth 
declines  are  very  similar  in  timing  to  studies 
done  in  the  Northeast,  but  they  are  not  as 
geographically  widespread  and  are  less  consistent 
within  a  given  sample.  Additionally,  the 
analysis  of  tree  ring  data  is  extremely 
complicated  because  of  the  effects  of  tree  age, 
stand  competition,  climate,  and  physiological 
responses  to  stress  that  may  persist  for  several 
years.  Therefore  the  interpretation  of  these 
data  has  been  the  subject  of  considerable 
controversy. 

Vast  differences  exist  in  species'  composition, 
structure,  and  stand  productivity  within  the 
limited  extent  of  the  southern  spruce -fir.  These 
differences  are  a  result  of  climatic  disparities 
associated  with  elevational  gradients,  past 
management  histories,  and  other  site-specific 
variables.  These  environmental  factors,  and  the 
lack  of  historical  data  in  the  South,  preclude 
any  analysis  designed  to  discover  either  the 
causes  of  observed  declines  or  even  if  the 
observed  declines  are  abnormal  in  these  forests. 
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A  Tree  Ring  Analysis  of  Red  Spruce  in  the  Southern  Appalachian  {fountains 

Edward  R.  Cook 


INTRODUCTION 

A  recent  analysis  of  red  spruce  (Picea  rubens 
Sarg.)  tree  ring  widths  in  southern  Appalachian 
stands  has  caused  concern  that  the  red  spruce 
forests  of  the  southern  Appalachian  mountains  may 
be  in  an  early  stage  of  decline  (Adams  and  others 
1985).  Although  many  hypotheses  have  been 
generated  regarding  the  cause  of  the  red  spruce 
decline,  no  definite  answer  has  yet  been  found 
(McLaughlin  1985).  If  the  decline  in  red  spruce 
ring  width  can  be  explained  by  natural  effects , 
then  costly  and  needless  pollution  controls  may 
be  avoided. 

The  Tennessee  Valley  Authority  (TVA)  and  the 
National  Park  Service  (NPS)  conducted  studies  on 
permanent  plots  established  throughout  forests  of 
the  southern  Appalachians.  Long-term  changes  in 
the  composition  and  health  of  the  forest  were 
monitored.  Using  the  data  provided  by  the  NPS 
and  TVA,  the  objective  of  my  analysis  was  to 
determine  if  the  recent  patterns  in  the  ring 
widths  indicate  an  anomalous  decline  and  if  this 
decline  can  be  explained  by  natural  environmental 
factors  related  to  climate.  The  analyses  were 
performed  on  annual  tree  ring  chronologies 
(Fritts  1976,  Cook  1985)  developed  from  the  ring 
width  series  of  each  plot. 

TREE  RING  DATA  QUALITY  CHECK  AND  STANDARDIZATION 

The  quality  of  the  data  provided  by  the  NPS  and 
TVA  was  checked  using  the  COFECHA  program  of 
Holmes  (1983).  The  program  checked  for  cross- 
dating  errors,  measurement  errors,  and  other  ring 
width  irregularities  that  might  limit  the 
usability  of  ring  width  time  series  for  tree  ring 
analysis.  Because  the  actual  increment  cores 
were  not  available  for  this  quality  check,  the 
program  output  was  used  to  verify  cross -dating, 
make  corrections  of  dating  when  possible,  and 
eliminate  ring  width  series  for  which  no  obvious 
corrections  could  be  made.  Approximately  15 
percent  of  the  ring  width  series  were  either 


corrected  or  removed  from  the  data  set. 
Therefore  the  number  of  ring  width  series  for 
some  plots  were  reduced  to  as  few  as  three. 

After  the  quality  check,  the  remaining  ring 
width  series  of  each  plot  were  standardized 
(Fritts  1976,  Cook  1985)  to  remove  long-term 
trends  in  growth  associated  with  tree  age,  size, 
and  stand  dynamics.  The  need  for  standardization 
prior  to  creating  a  stand- average  tree  ring 
chronology  is  discussed  in  detail  in  Fritts 
(1976)  and  Cook  (1987).  Because  the  ring  width 
series  were  rarely  more  than  100  years  long, 
negative  exponential  or  linear  regression  curves 
were  used  to  detrend  the  series . 

However,  it  was  unlikely  that  this  conservative 
de trending  method  would  remove  any  anomalous 
decline  signal  during  standardization.  After  the 
growth  curve  was  estimated  for  each  series,  the 
tree  ring  indices  were  computed  as: 

It  =  Rt/  Gt 

where  It  equaled  the  tree  ring  index,  Rt  was  the 
ring  width,  and  Gt  equaled  the  growth  curve 
value,  all  for  year  t.  Therefore  a  tree  ring 
index  can  be  defined  as  the  ratio  of  the  actual 
ring  width  to  the  expected  value  as  estimated  by 
Gt.  Tree  ring  indices  have  a  long-term  mean  of 
1.0  and  a  variance  that  is  reasonably  time 
stable.  Thus  tree  ring  indices  are  stationary 
processes  that  can  be  averaged  into  a  stand- 
average  series.  After  each  ring  width  series  was 
reduced  to  index  form,  the  tree  ring  index  series 
of  each  plot  were  averaged  into  a  final  tree  ring 
chronology  using  the  biweight  robust  mean 
(Mosteller  and  Tukey  1977)  to  reduce  the 
influence  of  outliers  on  the  computation  of  the 
mean- value  function. 

STRATIFICATION  AND  SCREENING  OF  PLOTS 

The  TVA  and  NPS  plots  were  stratified  by 
elevation  into  three  groups:  below  5,400  feet, 
5,400  to  6,000  feet,  and  above  6,000  feet.  These 
strata  reflected  a  vegetational  gradient  in  the 
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mountains.  The  lowest  red  spruce  stratum  was 
predominantly  a  mixed  conifer -hardwoods  forest 
zone,  while  the  strata  above  5,400  feet  were 
within  the  spruce -fir  zone  (White  1984) .  The 
highest  elevational  stratum  contained  the  highest 
proportion  of  Fraser  fir  (Abies  fraseri  [Pursh] 
Poir.)  relative  to  red  spruce  and  was  the  zone 
most  heavily  impacted  by  the  balsam  woolly 
adelgid  (Eagar  1984) . 

The  elevational  strata  also  reflected  a 
climatic  gradient  of  increasing  precipitation  and 
decreasing  temperature  with  increasing  elevation. 
Thus  the  below- 5 ,400- foot  stratum  was  the  warmest- 
and  driest  area  and  the  above -6 ,000- foot  stratum 
the  coolest  and  wettest.  Additionally,  the 
highest  stratum  was  enveloped  in  clouds  most 
often  and  the  lowest  stratum  was  enveloped  in 
clouds  least  often. 

The  climatic  response  of  red  spruce  in  the 
southern  Appalachian  Mountains  has  not  been 
studied  as  thoroughly  as  it  has  in  the  northern 
Appalachians  (Conkey  1979,  McLaughlin  and  others 
1987,  Cook  and  others  1987).  However,  based  on 
ecological  principles,  it  is  probable  that  the 
role  of  temperature  will  increase  as  a  limiting 
factor  in  red  spruce  growth  at  the  elevational 
extremes  of  the  species'  range.  Precipitation 
will  also  be  more  important  as  a  limiting  factor 
at  the  below- 5, 400 -foot  plots.  Therefore  a 
gradient  in  the  response  of  red  spruce  to  climate 
should  be  found  that  will  correlate  well  with 
some  aspect  of  the  known  climatic  gradient. 

The  stratification  by  elevation  produced  11  TVA 
and  8  NPS  plots  in  the  above -6 , 000-foot  stratum, 
13  TVA  and  6  NPS  plots  in  the  5,400-  to  6,000- 
foot  stratum,  and  24  TVA  and  7  NPS  plots  in  the 
below- 5 ,400 -foot  stratum.  At  this  stage,  some  of 
the  plots  were  eliminated  because  of  the 
shortness  of  the  tree  ring  chronologies.  Because 
1930  was  chosen  as  an  initial  criterion  for  the 
inclusion  of  tree  ring  series  in  the  analyses, 
any  series  beginning  after  1930  was  eliminated. 
The  year  1930  allowed  for  the  inclusion  of  the 
large  majority  of  plots  and  simultaneously 
provided  an  adequate  time  base  for  the 
dendroclimatic  analyses.  A  longer  time  base 
would  have  been  better,  but  it  also  would  have 
eliminated  too  many  sites.  Furthermore,  the  best 
available  climatic  data  began  in  1931. 

After  the  elimination  of  short  series,  the 
sample  depth  of  the  1930  decade  was  examined  for 
each  remaining  chronology.  If  the  sample  depth 
was  largely  based  on  only  one  increment  core  in 
that  decade  ,  then  that  chronology  was 
eliminated.  The  reason  for  this  step  in  the 
preliminary  screening  was  to  insure  that  the  1930 
decade  would  not  be  unduly  affected  by  poor 
replication  at  some  sites. 

The  final  result  of  the  screening  was  the 
selection  of  13  plots  above  6,000  feet,  15  plots 
between  5,400  and  6,000  feet,  and  29  plots  below 
5,400  feet.  Eight  above-6 ,000 -foot  series  were 
from  NPS  plots,  four  were  from  TVA  Mt .  Mitchell, 
and  one  was  from  TVA  Roan  Mountain.  Thus  the 
above-6 ,000 -foot  series  comprised  68  percent  of 
the  available  plots.  The  geographic  coverage  of 
this  stratum  was  obviously  limited  by  the  maximum 
elevations  of  the  mountains.  The  5,400-  to 
6, 000- foot  stratum  was  composed  of  six  NPS  plots 
and  nine  TVA  plots;  these  represented  78  percent 
of  the  available  plots.   The  geographic  coverage 


was  much  better.  Only  Roan  Mountain  and 
Grandfather  Mountain  were  not  represented.  The 
below- 5 ,400- foot  stratum  was  composed  of  7  NPS 
and  22  TVA  plots,  representing  93  percent  of  the 
available  plots.  The  geographic  coverage  was 
complete,  with  all  mountains  represented. 

AUTOREGRESSIVE  TIME  SERIES  MODELING 

Tree  ring  series  invariably  possess  some  degree 
of  serial  persistence  or  autocorrelation  that  is 
principally  due  to  physiological  preconditioning 
within  the  tree .  Therefore  the  information 
contained  in  a  given  ring  width  is  somewhat 
determined  by  past  tree  growth  and  vigor. 
Typically,  the  autocorrelation  structure  of  tree 
ring  indices  can  be  adequately  modeled  as  an 
autoregressive  process  (Cook  1985) .  The  general 
autoregressive  (AR)  process  of  order  p  has  the 
form  (Box  and  Jenkins  1970) : 
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where  Zt  is  the  observed  process  for  year  t,  et 
equals  an  unobserved  input  or  random  shock  that 
does  not  contain  any  autocorrelation,  and  *^  the 
autoregressive  coefficients  of  the  AR(p)  process. 

In  the  context  of  this  tree  ring  analysis,  the 
Zt  were  the  tree  ring  indices  for  a  plot.  Each 
tree  ring  series  was  modeled  as  an  AR(2)  for  the 
common  interval  1930-83.  The  choice  of  an  AR(2) 
model  was  based  on  previous  experience  modeling 
longer  red  spruce  chronologies  as  AR  processes. 
The  common  AR  persistence  structure  among  all 
plots  within  each  stratum  was  also  estimated 
using  a  pooling  procedure  described  in  Cook 
(1985) .  Differences  between  the  common  AR  model 
and  those  for  the  individual  series  were  useful 
tools  for  measuring  the  level  of  autocorrelated 
noise  in  the  individual  series,  which  may  have 
been  caused  by  different  stand  histories  and 
disturbance  regimes. 

For  the  above-6 ,000- foot  stratum,  the  common  or 
pooled  AR  coefficients  and  the  percent  variance 
explained  by  autoregression  (R  )  were:  $^  = 
0.461,  *2  =  0.284,  and  R2  =  46.1  percent.  For 
the  13  individual  series,  the  average  statistics 
were:  *i  =  0.566,  $2  =  0.138,  and  R2  =  47.9 
percent.  Although  the  R2's  were  similar,  the  AR 
coefficients  were  noticeably  different,  probably 
because  of  residual  trend  or  trendline  lack-of- 
fit  in  the  individual  series.  However,  the 
similarity  of  the  R2 '  s  suggested  that  the 
differences  between  the  tree  ring  chronologies 
were  largely  random  through  time.  Thus  the  long- 
term  disturbance  histories  of  these  plots  may 
have  been  similar  since  1930. 

For  the  5,400-  to  6,000-foot  stratum,  the 
pooled  statistics  were:  *^  =  0.324,  $2  =  0.177, 
and  R2  -  18.1  percent.  For  the  15  individual 
series,  the  average  statistics  are:  §\  =  0.484, 
*2  ~  0.168,  and  R2  =  39.1  percent.  The  pooled 
AR(1)  coefficient  and  R2  were  considerably 
smaller  than  the  average  values  for  the 
individual  series.  The  latter  indicated  a  high 
level  of  autocorrelated  noise  or  out-of-phase 
behavior  between  series.  Therefore  the 
disturbance  histories  of  these  plots  were 
probably  more  variable  than  those  in  the  higher 
stratum. 

For  the  below- 5 ,400- foot  stratum,  the  pooled 


statistics  were:  *i  =  0.422,  $2  "  0.124,  and  R2 
-  24.4  percent.  For  the  29  individual  series, 
the  average  statistics  were:   $^  =  0.530,  $2 


0.094,  and  R 


2  . 


39.3  percent.   These  statistics 


were  close  to  those  from  the  intermediate 
elevation  plots  and  indicated  a  similar  degree  of 
nonhomogeneity  from  plot  to  plot. 

The  lower  levels  of  plot  homogeneity  in  the 
strata  below  6,000  feet  suggested  that  these 
plots  have  more  varied  disturbance  histories.  An 
examination  of  the  individual  time  series  from 
these  plots  confirmed  this  inference.  Some  of 
the  plots  showed  release  patterns  early  in  this 
century  that  were  consistent  with  logging 
activity.  Given  the  much  reduced  spatial 
coverage  of  the  above -6 ,000 -foot  plots,  the 
higher  level  of  homogeneity  among  these  plots  was 
probably  related  to  the  lack  of  interference  by 


PRINCIPAL  COMPONENTS  ANALYSIS  (PCA) 

Because  each  chronology  was  based  on  a  small 
sample  of  trees,  the  dendroclimatic  modeling  of 
"  ~h  plot  chronology  was  not  considered  a  viable 
approach.  The  results  would  have  been  somewhat 
chaotic  because  of  the  very  high  level  of  noise 
in  each  chronology.  Therefore  the  common 
variance  among  all  series  within  each  stratum  was 
pooled  using  principal  components  analysis  (PCA) 
(Cooley  and  Lohnes  1971).  In  PCA,  the  structure 
in  the  correlation  matrix  of  variables  is 
transformed  into  a  new  set  of  uncorrelated  or 
orthogonal  modes  of  behavior  called 
eigenvectors.  Each  eigenvector  accounts  for  a 
unique  proportion  of  the  cotal  variance  in  the 
original  data.  The  first  eigenvector  associated 
with  the  largest  eigenvalue  accounts  for  the 
greatest  percentage  of  common  variance  among  all 
variables  in  the  correlation  matrix. 

Each  eigenvector  is  composed  of  a  number  of 
loadings  or  coefficients  equal  to  the  number  of 
original  variables  in  the  correlation  matrix. 
These  loadings,  which  may  be  positive  or 
negative,  reflect  the  relationships  between 
variables  for  a  specific  eigenvector.  Frequently 
the  loadings  of  the  first  eigenvector  are  all 
positive  or  negative,  which  indicates  that  the 
variables  being  analyzed  all  behave  similarly. 
Therefore,  in  this  study,  the  tree  ring  series  at 
each  elevational  stratum  could  exhibit  a  common 
signal  due  to  climate,  disturbance,  or  pollution. 

The  loadings  of  the  first  eigenvector  can  also 
be  used  to  create  a  time  series  of  scores  that 
reveal  how  this  most  common  component  among  all 
series  behaves  through  time.  This  series  of 
scores  is  similar  to  a  weighted  mean,  because 
each  series  is  weighted  by  its  eigenvector 
loading  and  then  summed  with  the  other  weighted 
series  for  each  year.  The  weighting  scheme  is 
optimal  in  the  sense  that  no  other  component  can 
account  for  more  of  the  common  variance  between 
series  than  the  first  eigenvector.  Consequently, 
the  scores  for  each  elevational  stratum  should 
have  a  strong  common  signal  for  dendroclimatic 
analysis . 

The  PCA  analyses  were  done  twice  for  each 
stratum:  once  on  the  original  tree  ring  indices 
and  again  on  the  indices  after  removing  AR(2) 
persistence  from  each  series.  Indices  are 
referred  to  as  prewhitened  after  removal  of  AR(2) 
persistence.   For  the  above-6 ,000 -foot  stratum  of 
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13  plots,  the  first  eigenvector  of  the  original 
tree  ring  indices  accounted  for  50.7  percent  of 
the  total  variance,  while  that  of  the  prewhitened 
indices  accounted  for  53.6  percent  of  the 
variance .  Common  variance  increased  after 
prewhitening  because  of  a  reduction  of  noise 
variance  resultant  from  autoregressive  modeling. 
In  figure  1,  the  loadings  for  this  eigenvector 
are  all  positive,  indicating  an  existing  common 
signal  among  all  series.  The  loadings  for  the 
prewhitened  indices  were  more  uniformly  positive 
than  those  of  the  original  indices.  This 
uniformity  indicated  that  some  of  the  differences 
between  the  original  indices  were  amplified  by 
the  autoregression  within  those  series. 

For  the  5,400-  to  6, 000 -foot  stratum  of  15 
plots ,  the  first  eigenvector  of  the  original  tree 
ring  indices  accounted  for  34.1  percent  of  the 
variance,  while  the  prewhitened  indices  accounted 
for  45.9  percent.  The  larger  increase  in  common 
variance  after  prewhitening  indicated  that  these 
chronologies  were  less  homogeneous  than  those  in 
the  higher  stratum.  Greater  variability  in  site 
characteristics  and  stand  histories  in  this 
intermediate  stratum  may  have  caused  the 
difference.  Comparing  the  single  series  and 
pooled  autoregression  models  also  indicated  the 
lack  of  homogeneity  between  the  chronologies. 
However,  the  more  restricted  geographic  coverage 
of  the  high  stratum  may  be  a  biasing  agent  in 
this  comparison.  As  before,  the  eigenvector 
loadings  (fig.  1)  were  also  more  uniform  after 
prewhitening. 

For  the  below- 5 ,400- foot  stratum  of  29  plots, 
the  first  eigenvector  of  the  original  tree  ring 
indices  accounted  for  32.6  percent  of  the 
variance ,  while  that  of  the  prewhitened  indices 
accounted  for  40.7  percent.  The  magnitude  of  the 
difference  was  similar  to  that  of  the 
intermediate  stratum.  Therefore  the  level  of 
homogeneity  between  plots  was  similar,  an 
inference  also  supported  by  the  earlier 
autoregressive  modeling  results.  The  eigenvector 
loadings  (fig.  1)  of  the  prewhitened  indices  were 
also  more  uniform  than  the  loadings  of  the 
original  indices. 

Generally  the  strength  of  the  common  signal 
within  each  stratum  was  directly  correlated  with 
the  elevational  gradient.  The  tree  ring  patterns 
of  high  plots  were  more  similar  among  themselves 
than  those  of  the  lower  plots .  Although  it  may 
appear  that  this  result  reflected  more  limiting 
growth  conditions  towards  the  upper  elevational 

limit,  the  bias  in  the  geographic  coverage  of 
that  stratum  limits  the  strength  of  this 
interpretation. 

The  eigenvector  amplitudes  or  scores  of  each 
stratum  are  shown  in  figure  2.  The  solid  line 
plots  were  derived  from  the  orignal  tree  ring 
indices ,  and  the  dashed  line  plots  were  derived 
from  the  AR(2)  prewhitened  indices. 

The  scores  derived  from  the  original  indices 
indicate  an  overall  pattern  of  below- average 
growth  at  all  plots  since  about  1966.  The 
largest  departure  was  for  the  above -6, 000 -foot 
plots.  The  average  score  since  1966  was  -2.65 
with  a  standard  error  of  ±0.495.  For  the  5,400- 
to  6,000-foot  stratum,  the  average  score  was- 
1.54  ±0.522.  And,  for  the  below-5 , 400-foot 
stratum,  the  average  score  was  -1.73  ±0.638. 
These  long-term  departures  appeared  to  exceed  the 
95-percent  significance  level  using  a  simple  t- 
test.   However,   the   use   of   a   t-test   on 
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Figure  1 . - -The  eigenvector  loadings  of  the  first  principal  component  of  each 
red  spruce  elevation  stratum.  The  solid  line  plots  correspond  to 
the  loadings  of  the  original  tree  ring  index  chronologies .  The 
dashed  line  plots  correspond  to  the  loadings  of  the  same  series 
after  second-order  autoregression  has  been  removed  from  each  one. 
The  eigenvectors  were  extracted  from  the  correlation  matrix.  The 
percent  variance  accounted  for  by  each  eigenvector,  original 
(prewhitened) ,    is    indicated . 


ABOVE  8000-FT  BLEV.  PLOTS 1ST  AMPLITUDE    50.7X  (53.6%)  VARIANCB 


5400-6000 -FT  ELEY.  PLOTS 1ST  AMPLITUDE    34. IX  (45.9%)  VARIANCE 


BELOW-5400-PT  ELEV.  PLOTS 1ST  AMPLITUDE    32.6%  (40.7%)  VARIANCE 
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Figure  2. --The  principal  component  amplitudes  or  scores  corresponding  to  the 
eigenvectors  in  figure  1.  The  solid  line  plots  are  for  the 
original  series.  The  dashed  line  plots  are  for  the  prewhitened 
series . 
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autocorrelated  time  series  such  as  these  can  be 
extremely  misleading.  The  number  of  independent 
observations  and  the  degrees-of -freedom  may  be 
much  less  than  the  number  indicated  by  the 
available  observations. 

One  way  of  avoiding  the  negative  effect  of 
autocorrelation  on  the  degrees-of -freedom  is  to 
use  the  scores  of  the  AR(2)  prewhitened  indices. 
By  definition,  these  scores  do  not  have  any 
significant  autocorrelation  related  to  that  level 
of  autoregression.  In  figure  2,  the  prewhitened 
scores  indicate  a  smaller  reduction  in  growth 
since  1966  in  all  series.  The  post-1965  means 
confirmed  this.  For  the  above- 6 ,000 -foot ,  5,400- 
to  6,000-foot,  and  below- 5 ,400- foot  strata,  the 
1966  to  1983  means  were  -1.09  ±0.779,  -0.42 
±0.807,  and  -0.71  ±0.992,  respectively.  None  of 
these  means  passed  a  t-test  at  the  95 -percent 
significance  level.  Therefore  there  may  have  not 
been  any  reduction  in  growth  since  1966.  The 
decline  in  the  original  tree  ring  indices  may  be 
largely  explained  by  the  endogenous 
autoregressive  persistence  of  the  tree  ring  data 
and  the  way  in  which  it  amplifies  the  behavior  of 


the   random   shocks , 


which   are   largely 


exogenous  to  the  plots. 

The  above  conclusion  was  conservative  because 
the  AR  coefficients  used  for  prewhitening  were 
based  on  information  in  both  the  pre -1966  and 
post- 1966  time  periods.  The  method  employed 
minimized  the  probability  of  a  type  I  error 
because  it  assumed  that  the  AR  coefficients  had 
not  changed  through  time  despite  an  intervention 
in  the  et  that  may  have  occurred  in  1966.  Since 
an  intervention  in  the  et  could  have  a  strong 
impact  on  the  estimation  of  the  AR  coefficients, 
the  prewhitening  may  have  removed  part  of  the 
response  to  an  intervention  had  it  occurred.  In 
order  to  reduce  the  probability  of  a  type  II 
error  in  these  analyses,  an  alternate  method  of 
intervention  analysis  (Box  and  Tiao  1975)  testing 
for  the  occurrence  of  an  intervention  in 
autocorrelated  time  series  was  used. 

INTERVENTION  ANALYSIS 

Intervention  analysis  specifically  allows  for 
autocorrelation  when  testing  for  the  occurence  of 
an  intervention  in  time  series.  A  simple  form  of 
intervention  analysis  was  used  in  this  study  to 
test  for  the  occurrence  of  a  decline  in  the 
scores  since  1966.  The  form  of  the  intervention 
chosen  was  a  step -function,  which  is  expressed  as 
[0  0  0  ...]  from  1930  to  1965  and  [1  1  1  ...] 
from  1966  to  1983.  This  step-function  served  as 
one  of  the  predictor  variables  in  the  analysis. 
To  account  for  autocorrelation  in  each  time 
series,  the  scores  for  years  t-1  and  t-2  were 
also  used  as  predictors.  The  model  used  allowed 
for  both  the  occurrence  of  a  step  reduction  in 
growth  and  AR(2)  persistence.  The  intervention 
model  was  set  up  as  a  multiple  regression 
analysis  problem.  In  each  case,  only  the  step- 
function  and  the  lag-1  variable  proved  to  be 
statistically  significant  at  the  90-percent  level 
or  higher. 

In  contrast,  the  lag- 2  variable  never  exceeded 
the  60-percent  significance  level.  For  this 
reason,  lag- 2  was  not  used  in  the  final  models. 
Although  one  might  infer  that  the  previous  AR(2) 
models   were   reasonable   only   because   an 


intervention  around  1966  changed  the  system,  the 
timing  of  the  intervention  was  hypothesized  only 
after  an  examination  of  the  data.  Thus  any 
inferences  concerning  a  change  in  persistence 
structure  because  of  an  intervention  must  take 
into  account  the  a  posteriori,  nature  of  these 
analyses.  This  issue  will  be  addressed  later,  as 
it  affects  significance  tests. 

The  results  of  the  intervention  analysis  were 
as  follows : 


Stratum 


ARm 


Step 


R 


T 


Above  6,000  ft    0.260 
5,400-6,000  ft    0.234'' 
Below  5.400  ft    0.360'' 


** 


-0.566' 
-0.369* 
-0.275 


** 


58.9% 
27.6% 
28.8% 


p  <0.10, 


**p  <0.05,      ***  p  <0.01 


The  strength  of  the  step  intervention  was 
directly  correlated  with  elevation.  The  above- 
6, 000- foot  scores  showed  the  strongest  indication 
of  an  intervention  in  1966,  which  resulted  in  a 
steplike  reduction  in  growth.  This  result  was 
consistent  with  the  original  examination  of  the 
1966  to  1983  means  for  these  scores.  However, 
the  step -elevation  relationship  was  new.  The 
high  level  of  persistence  in  the  above -6 ,000- foot 
scores  was  greatly  reduced  by  the  step.  In 
contrast,  the  persistence  in  the  lower  strata  was 
reduced  less.  The  reduction  in  persistence  from 
the  earlier  AR(2)  modeling  appears  to  be 
proportional  to  the  strength  of  the  step. 

The  probability  levels  of  the  intervention 
analysis  were  based  on  an  a  priori  significance 
test  in  each  case.  Acceptance  of  these  results 
would  effectively  minimize  type  II  error  at  the 
expense  of  type  I  error,  in  contrast  to  the 
earlier  prewhitening  results  that  minimized  type 
I  error.  Therefore  these  sets  of  results  served 
as  useful  limits.  As  noted  earlier,  there  was  a 
problem  in  applying  a  priori  significance  tests 
to  a  statistical  analysis  problem  that  was 
principally  based  on  an  a  posteriori  examination 
of  the  data.  Furthermore,  the  a  posteriori 
examination  of  the  scores  for  an  intervention 
allowed  for  50  possible  intervention  dates  for  a 
step -function. 

Based  on  probability  theory,  the  probability  of 
finding  a  statistically  significant  step- 
function  under  such  conditions  is  related  to  the 
a  priori  significance  level  as: 

P  =  1  -  (1  -  p)m 

where  P  is  the  a  posteriori  probability  level,  p 
is  the  a  priori  probability  level  of  the  test 
being  applied,  and  m  is  the  number  of  times  the 
test  could  be  applied  to  the  data.  If  this 
correction  is  applied  to  the  probability  levels 
for  the  step  interventions  shown  above,  only  the 
above  -  6000- foot  step  intervention  remains 
statistically  significant  (P<0.01).  The  other 
steps  do  not  even  pass  the  P<0.50  level.  This 
correction  is  probably  overly  severe  since  the  a 
priori  information  about  the  probable  timing  of 
red  spruce  decline  in  the  northern  Appalachians 
(Johnson  and  Siccama  1983) . 

All  available  evidence  indicated  that  the 
decline  in  the  northern  Appalachians  started  in 
the  late  1950'  s  to  early  1960's.   There  is  no 
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evidence  to  suggest  that  any  decline  in  the 
southern  Appalachians  began  before  the  northern 
decline.  Thus  no  intervention  should  be  found 
prior  to  1960.  If  this  information  is  used  to 
limit  the  intervention  time  window  to  1960-81, 
the  results  of  the  lower  two  strata  still  remain 
well  outside  the  P<0.10  level  of  significance, 
which  is  still  unacceptable.  For  an  a  posteriori 
probability  level  of  P=0.10  to  be  achieved,  an  a 
priori  probability  of  p=0.01  and  an  intervention 
time  window  of  10  years  are  needed.  Unless 
additional   constraints   based   on   a   priori 


information  can  be  found  to  reduce  the  time 
window  of  the  hypothesized  intervention,  the  null 
hypothesis  of  no  intervention  for  the  lower 
strata  cannot  be  rejected  on  statistical  grounds. 

Figure  3  shows  each  series  of  scores  with  its 
fitted  intervention  model.  As  the  statistics 
indicate,  the  above -6 ,000 -foot  model  revealed  a 
much  more  pronounced  step  change  than  the  other 
models . 

All  strata  showed  signs  of  decreasing  growth 
after  1966.  However,  growth  reduction  in  the 
lower   strata   diminished   with   decreasing 
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Figure  3. --The  actual  (solid)  and  estimated  (dash)  tree  ring  scores  based  on 
fitting  a  step  intervention  model  to  each  series.  The  date  of 
the   intervention   is   1966. 
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elevation.  The  cause  of  the  growth  reduction  is 
still  undetermined.  Climate  or  other  changes  in 
natural  influences  are  as  probable  a  cause  as  are 
anthropogenic  pollutants. 

The  reality  of  the  elevational  gradient  in  step 
size  may  be  questioned.  However,  the  existence 
of  known  environmental-climatic  gradients  in  the 
mountains  suggests  that  the  step-size  gradient  is 
a  reality.  The  step- size  gradient  may  be  due  to 
a  temperature- related  phenomenon.  This 
hypothesis  is  consistent  with  what  is  now  known 
about  the  relationship  between  temperature  stress 
and  red  spruce  declines  in  the  northern 
Appalachian  Mountains  (Cook  and  others  1987). 


DENDROCLIMATOLOGY  OF  RED  SPRUCE 
BY  ELEVATIONAL  STRATUM 

As  noted  earlier,  the  three  strata  used  in  this 
study  follow  both  vegetational  and  climatic 
gradients,  which  are  directly  correlated  with 
elevation.  To  illustrate  the  reflection  of  this 
gradient  in  the  tree  rings ,  figure  4  shows  the 
three  series  of  original  tree  ring  index  scores 
superimposed  on  each  other.  Various  time  periods 
in  figure  4  (1932-34,  1935-57,  1959-63)  indicate 
striking  gradients  across  scores  correlated  with 
elevation.  The  presence  of  these  gradients 
across  scores  suggests  that  an  elevational 
gradient  in  the  climatic  response  of  red  spruce 
operates  at  times . 

At  other  times  in  the  scores  (1940-42,  1958, 
1969) ,  the  gradient  breaks  down,  and  the  scores 
of  all  strata  are  similiar.  This  similarity 
suggests  that  the  relationship  between  climatic 
response  and  elevation  is  nonstationary  through 


time.  The  degree  to  which  the  gradient  exists 
probably  depends  on  which  climatic  variables  are 
limiting  to  red  spruce  growth  in  a  given  year  and 
how  those  climatic  variables  are  influenced  by 
elevation.  For  example,  based  on  the  physics  of 
precipitation  formation  and  its  interaction  with 
orography,  the  influence  of  drought  on  red  spruce 
growth  should  diminish  with  increasing  elevation. 
However,  once  the  available  moisture  supply  is  no 
longer  limiting  to  growth,  this  drought  response 
gradient  would  probably  disappear  from  the  tree 
rings . 

In  this  study,  dendroclimatic  modeling  is 
limited  by  the  lengths  of  the  series  being 
modeled  and  the  unavailability  of  long  climatic 
time  series.  Ideally,  the  modeling  should 
proceed  as  described  by  Cook  1987;  the 
dendroclimatic  signal  should  be  modeled  for  a 
long  preintervention  time  period  of  perhaps  50  to 
60  years.  A  model  should  then  be  used  to 
forecast  or  predict  tree  rings  through  a  period 
to  the  present  that  includes  both  another 
preintervention  time  block  and  the  post- 
intervention  period.  The  time  stability  of  the 
dendroclimatic  model  must  be  tested;  therefore 
another  preintervention  time  period  is  needed. 
If  the  model  is  verified  as  time  stable,  then  it 
can  be  used  to  test  for  an  intervention  that 
changes  the  tree  ring  response  to  the  model. 
This  method  has  been  successively  used  in 
analyzing  the  red  spruce  decline  in  the 
Appalachian  Mountains  (McLaughlin  and  others 
1987,  Cook  and  others  1987). 

Because  of  the  insufficient  tree  ring  time 
base,  a  weaker  method  of  modeling  was  implemented 
that  provided  a  basis  for  inference  regarding  the 
climatic  response  gradient  hypothesized  earlier. 
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Figure  4 . - -The  amplitudes  of  the  original  tree  ring  indices  superimposed  on 
each  other.  The  purpose  of  this  plot  is  to  highlight  certain 
time  periods  when  elevation-related  gradients  in  climate  response 
are   likely   to   be   occurring. 
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This  method  was  based  on  simple  correlations 
between  the  tree  ring  scores  and  monthly  climatic 
data  for  the  period  1931-83. 

Prior  to  the  correlation  analyses,  the  tree 
ring  scores  for  each  elevational  stratum  were 
prewhitened  to  remove  autocorrelation.  In  each 
case,  significant  AR(1)  or  AR(2)  persistence  was 
removed.  The  monthly  temperature  and 
precipitation  data,  averaged  over  the  northern 
and  southern  mountain  climatic  divisions  of  North 
Carolina  and  the  southwestern  mountain  division 
of  Virginia,  were  similarly  modeled  for 
autocorrelation.  In  this  case,  the  climatic  data 
showed  very  weak  or  nonexistent  autocorrelation 
out  to  lag  3 .  Therefore ,  the  climatic  data  were 
not  prewhitened. 

The  dendroclimatic  modeling  was  then  treated  as 
a  multiple  input-single  output  transfer  function 
(Box  and  Jenkins  1970)  in  which  ring  width  was  a 
function  of  climate.  Given  the  lack  of 
autocorrelation  in  either  the  input  or  output 
series,  the  principal  aim  of  the  transfer 
function  model  was  to  identify  those  climatic 
variables  that  correlated  significantly  with  tree 
rings  and  identify  any  delay  or  lag- response 
between  the  inputs  and  the  output.  The  analysis 
assumed  that  the  climatic  variables  were 
orthognal,  an  assumption  that  was  violated  for 
almost  all  variables.  This  violation  may  have 
increased  the  number  of  significant  climatic 
variables  in  the  model.  However,  since  the  aim 
of  these  analyses  is  strictly  correlative  and  not 
predictive ,  this  should  not  have  any  serious 
impact  on  the  results . 

The  tree  ring  scores  were  lagged  up  to  3  years 
in  the  transfer  function  analyses,  meaning  that 
each  monthly  climatic  variable  was  correlated 
with  each  series  of  scores  for  years  t,  t+1,  t+2, 
and  t+3 .  A  plot  of  the  correlations  by  lag 
produced  a  normalized  form  of  the  impulse 
response  function  for  1932-62  and  1966-80  time 
periods  of  the  system  being  modeled  (Box  and 
Jenkins  1970).  For  each  period,  48  precipitation 
and  48  temperature  correlation  coefficients  were 
computed.  While  the  a  posteriori  nature  of  these 
analyses  makes  the  use  of  a  priori  significance 
tests  very  questionable,  these  results  should  be 
viewed  as  more  exploratory  than  confirmatory. 
Therefore  the  a  priori  confidence  limits  will  be 
used  to  assess  the  significance  of  the 
correlation  coefficients. 

The  results  of  this  modeling  were  somewhat 
complex  to  explain.  In  each  time  period,  some 
indications  of  climatic  gradients  were  found. 
For  example,  in  the  1932-62  period,  the 
correlation  between  tree  rings  and  March 
precipitation  at  lag  t+2  were: 

Above  6,000  feet:  -0.102 
5400  to  6,000  feet:  -0.604 
Below  5,400  feet:    -0.636 

The  correlations  of  the  lower  two  strata  were 
significant  (p<0.001)  in  a  statistical  sense. 
However,  the  t+2  lags  were  very  difficult  to 
explain  physiologically.  More  disconcerting, 
these  correlations  completely  lost  statistical 
significance  (maximum  | r |<0. 15)  in  the  1966-80 
period.  Therefore  these  significant  correlations 
were  either  spurious  or  the  climatic  signal  in 
the  red  spruce  was  highly  nonstationary.    The 


latter  problem  may  also  indicate  a  loss  of 
climatic  signal  comparable  to  what  has  apparently 
happened  to  the  declining  red  spruce  in  the 
northern  Appalachian  Mountains  (Cook  1987, 
McLaughlin  and  others  1987,  Cook  and  others, 
1987). 

In  the  suite  of  96  total  correlations ,  only 
three  monthly  temperature  variables  showed  any 
consistency  through  both  time  periods;  July, 
August,  and  September  temperatures  correlated 
with  t+1  lagged  tree  rings  as  follows : 


Stratum 

1932-62 

1966-80 

Above  6,000 

feet 

July 

-0.343* 

-0.273 

Augus  t 

-0.268 

-0.134 

Sept. 

-0.223 

-0.132 

5,400-6,000 

feet 

July 

-0.388** 

-0.395' 

August 

-0.334* 

-0.361 

Sept. 

-0.403** 

-0.326 

Below  5,400 

feet 

July 

-0.351* 

-0.5471 

August 

-0.247 

-0.576' 

Sept. 

-0.365** 

-0.453' 

*p<0.10**, 

p<0 . 05 

,   ***p<0 

.01 

There  is  an  indication,  especially  in  the  1966-80 
period,  of  an  elevational  gradient  in  the 
response  to  the  temperature  variables.  The  high- 
elevation  stands  seem  to  be  less  sensitive  to 
summer  temperatures  than  the  lower  stands .  There 
is  also  an  indication  that  the  below- 5, 400 -foot 
spruce  have  been  more  sensitive  to  summertime 
temperature  since  1966. 

On  the  basis  of  these  monthly  temperature 
correlations,  the  July,  August,  and  September 
temperatures  were  averaged  into  a  summer  season 
temperature  series  (fig.  5) .  Of  particular 
interest  is  the  summer  of  1980,  the  warmest 
summer  in  the  southern  Appalachians  since  1931. 
According  to  figure  4,  the  poorest  growth  year 
for  red  spruce  at  all  elevations  was  1981.  Given 
the  t+1  lag  response  of  red  spruce  to  summer 
temperatures  indicated  above,  the  poor  1981 
growth  year  was  probably  related  to  excessively 
warm  summer  temperatures  in  1980,  which  extended 
to  the  highest  elevations  in  the  mountains. 

Linear  regression  analyses  of  the  summer 
temperature  series  versus  prewhitened  red  spruce 
scores  indicated  that  the  below- 5 ,400-  and  5,400 
to  6, 000- foot  strata  were  equally  sensitive  to 
previous  summer  temperatures  over  the  period 
1932-83.  The  regression  R2's  were,  respectively, 
0.199  and  0.182.  In  contrast,  the  above -6,000 - 
foot  regression  R2  was  0.137.  The  prewhitened 
scores  and  their  temperature  estimates  are  shown 
in  figure  6.  Interestingly,  all  strata  follow 
the  pattern  of  summer  temperature  almost 
perfectly  since  1978.  This  corresponds  to  the 
warmer  than  average  temperatures  since  1977. 

It  was  previously  suggested  that  the  lower 
elevation  plots  should  be  more  stressed  by 
precipitation  deficiency  than  the  higher 
elevation  plots.  In  order  to  determine  the 
degree  of  drought  sensitivity  in  the  prewhitened 
red  spruce  scores ,  the  monthly  Palmer  Drought 
Severity  Indices  (PDSI)  (Palmer  1965)  were 
computed  from  the  divisional  average  temperature 
and  precipitation  data.   Simple  correlations  were 
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Figure  5. --July,  August,  and  September  average  temperatures  for  the  southern 
Appalachian  Mountains  since  1931.  The  dashed  line  highlights 
variance  at  frequencies  of  1/10  year  or  less. 


again  computed  between  monthly  PDSI  and  tree  ring 
scores  at  lags  t  through  t  +  3  for  the  time 
periods  1932-62  and  1966-80.  The  results  of 
these  simple  correlation  analyses  were  very 
similar  to  the  analyses  reported  earlier;  the 
correlations  were  generally  not  time  stable. 
However,  as  before,  summertime  drought  and  t  +  1 
lagged  scores  showed  some  time  stability  and 
statistical  significance. 

These  correlations  were  as  follows: 


Stratum 

Month 

1932-62 

1966-80 

Above   6 

000 

feet 

July 

0.066 

0.521** 

August 
Sept. 

-0.004 
-0.281 

0 . 640** 
0.728*** 

5,400-6 

000 

feet 

July 

August 

Sept. 

0.442** 
0.427** 
0.079 

0.521** 
0.666** 
0.761*** 

Below  5 

400 

feet 

July 

August 

Sept. 

0.457*** 

0.402** 

0.244 

0.458* 
0 . 641** 
0.771*** 

*p<0 . 10 ,   **p<0 . 05 ,   ***p<0 . 01 

In  the  1932-62  period,  there  was  a  clear 
indication  of  an  elevational  gradient.  The 
highest  stratum  showed  no  summertime  drought 
signal  in  contrast  to  the  lower  strata.  However, 
all  strata  showed  a  very  pronounced  summertime 
drought  response  in  the  1966-80  period.  This 
apparent  increase  in  sensitivity  to  PDSI  was 
stronger  than  that  indicated  by  summertime 
temperature  alone. 

As  before,  the  monthly  PDSI's  were  averaged 
into  a  summertime  season  estimate  of  drought 
since  1931  (fig.  7).   Of  particular  interest  is 


the  time  period  of  1952-55,  which  contained  the 
worst  drought  in  the  southern  Appalachians  since 
1931.  The  very  strong  elevation- related  gradient 
in  the  tree  ring  scores  for  this  period  was 
almost  definitely  caused  by  severity  of  this 
drought  and  the  way  in  which  it  diminished  with 
increasing  elevation.  It  is  difficult  to  explain 
why  red  spruce  at  all  elevations  showed 
approximately  the  same  level  of  response  to  PDSI 
since  1966.  Given  the  shortness  of  this  time 
period,  it  is  possible  that  these  results  are 
questionable,  even  with  the  high  significance 
levels  of  the  correlation  coefficients.  However, 
this  apparent  increase  in  sensitivity  to 
summertime  moisture  availability  should  be 
investigated  more  fully,  as  better  statistical 
methods  and  tree  ring  data  become  available. 

Linear  regression  analyses  of  the  prewhitened 
scores  versus  summer  PDSI  for  the  period  1932-83 
indicated  a  weaker  relationship  overall  than  for 


summer  temperature  alone . 


The  R^'s  for  the 
below-5,400-,  5,400-  to  6,000-,  and  above-6, 000- 
foot  strata  were  0.151,  0.114,  and  0.035, 
respectively.    The  actual  and  predicted  scores 

from  these  models  are  shown  in  figure  8.  There 
generally  appears  to  be  less  time  stability  in 
PDSI -spruce  relationships. 

Summer  temperatures  and  PDSI  were  correlated 
(r--0.38)  because  the  temperature  data  were 
partially  used  to  estimate  the  PDSI's.  However, 
the  level  of  correlation  was  not  high  enough  to 
indicate  that  the  PDSI  correlations  were 
completely  confounded  by  the  temperature  effects. 
In  fact,  when  the  PDSI  and  temperature  variables 
were  used  in  a  stepwise  multiple  regression 
analysis  to  predict  red  spruce  scores,  each 
variable  entered  the  model  according  to  the 
strength  and  sign  of  its  original  correlation 
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Figure   6. --Actual    (solid)   and  predicted    (dash)   prewhitened   tree  ring  scores. 
The   summer   temperature   series    (fig.    5)    was   used   as    the  predictor 
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Indices     (PDSI)  series    for     the    southern    Appalachian    Mountains. 

The  dashed  line  highlights  variance  at  frequencies   less   than  1/10 
year. 


with  the  scores.  The  resulting  R2 '  s  for  the 
below-5,400- ,  5,400-  to  6,000-,  and  above-6, 000- 
foot  strata  were  0.252,  0.203,  and  0.123, 
respectively.  The  actual  and  predicted  scores 
from  these  models  are  shown  in  figure  9 . 

SYNTHESIS  OF  RESULTS 

The  results  of  this  study  indicated  that  red 
spruce  in  the  southern  Appalachian  Mountains  have 
exhibited,  to  varying  degrees,  some  irregular 
behavior  in  their  ring  widths  since  the  mid- 
1960'  s.  At  elevations  above  6,000  feet, 
statistical  evidence  suggest  a  steplike  decline 
in  radial  increment  since  about  1966.  This 
decline  has  not  been  correlated  with  any  specific 
climatic  deviation  in  this  study.  However,  the 
way  in  which  the  magnitude  of  the  decline 
increased  with  elevation  suggested  that  the  cause 
of  the  decline  was  somewhat  related  to  elevation. 
A  more  thorough  search  for  natural  and 
anthropogenic  causes  of  this  putative  decline  is 
warranted.  In  addition,  new  and  improved 
collections  of  ring  width  data  are  highly 
desirable  to  refine  the  statistical  analyses  and 
validate  or  refute  the  intervention  results 
presented  here. 

The  decline  of  the  southern  red  spruce  at  high 
elevations  could  lead  to  broad  scale  mortality, 
as  found  in  northern  Appalachian  stands. 
However,  the  dendroclimatic  modeling  has  revealed 
an  apparent  singular  difference  between  the 
northern  red  spruce  and  southern  red  spruce 
conditions  in  the  Appalachian  Mountains  since  the 
1960's.  In  the  northern  red  spruce,  the 
dendroclimatic  signal  completely  disappeared 
after  the  trees  entered  the  post-1960  period  of 
declining  ring  widths  (Cook  1987,  McLaughlin  and 


others  1987,  Cook  and  others  1987).  However, 
based  on  the  temperature  modeling  demonstrated  in 
this  study  (fig.  6)  ,  the  dendroclimatic  signal 
appeared  to  continue  through  the  post- 1965 
decline  period.  Thus  the  reported  decline  does 
not  seem  to  represent  a  major  loss  of  tree 
vitality  as  was  indicated  for  the  northern  red 
spruce.  At  this  stage  of  inquiry,  the 
Appalachian  Mountain  northern  and  southern  red 
spruce  situations  appear  to  be  different. 

The  dendroclimatic  analyses  revealed  that 
previous  summer  temperatures  correlated 
significantly  with  red  spruce  ring  width  the 
following  year.  The  lag-1  negative  temperature 
correlations  were  remarkably  consistent  with  more 
rigorously  developed  dendroclimatic  models  for 
numerous  stands  of  red  spruce  in  the  northern 
Appalachians  (McGlaughlin  and  others,  [in  press], 
Cook  and  others  [in  press]).  It  is  increasingly 
clear  that  the  role  of  previous  summer 
temperature  as  a  determinant  of  red  spruce  growth 
and  vigor  is  genetically  based.  More 
importantly,  warm  summer  temperatures  appear  to 
be  strongly  correlated  with  past  and  present 
declines  of  red  spruce  in  the  northern 
Appalachians  (Cook  and  others  1987) .  Should  the 
apparent  increase  in  sensitivity  to  prior-summer 
temperatures  be  correct  for  the  below-5 , 400-foot 
spruce,  low-elevation  spruce  in  the  southern 
Appalachians  are  likely  to  decline  if  warmer  than 
average  summer  temperatures  persist.  A  warmer 
world  caused  by  CO2  and  other  greenhouse  gases 
would  not  positively  affect  the  future  of  red 
spruce  in  North  America. 

The  apparent  increase  in  red  spruce 
sensitivity  at  all  elevations  to  drought  since 
19§6  and  with  sensitivity  to  summer  temperatures 
since  1977,  suggests  that  the  southern  red  spruce 
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Figure  8 . - -Actual  (solid)  and  predicted  (dash)  prewhitened  tree  ring  scores. 
The  summer  PDSI  series  (fig.  7)  was  used  as  the  predictor  of  tree 
rings.      The  S?  of  each  model   is   indicated  by   the  RSQ  value. 
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■are  in  a  prolonged  period  of  climatic  stress.  A 
similar  pattern  of  increased  climatic  stress  from 
about  1938-60  preceeded  the  current  broad  scale 
decline  of  red  spruce  in  the  northern 
Appalachians  (Cook  and  others  1987).  Presently 
it  is  impossible  to  say  that  this  circumstantial 
agreement  in  symptomology  is  part  of  the 
epidemiology  of  red  spruce  decline.  However,  it 
is  cause  for  concern  and  warrants  further  study. 
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Utilizing  Time  Series  Models  and  Spatial  Analysis 
of  Forecast  Residuals  for  Tree  Ring  Analysis  of  Red  Spruce 


J.  Keith  Ord  and  Janice  A.  Derr 


SUMMARY 


METHODOLOGY  AND  RESULTS 


The  information  from  a  field  study  on  permanent 
plots  established  by  the  Tennessee  Valley 
Authority  in  the  Great  Smoky  Mountains  was  used 
to  detect  and  evaluate  recent  changes  in  annual 
ring  width  of  red  spruce  (Picea  rubens  Sarg.). 
Time  series  models  were  fit  to  mean  annual  ring 
widths  of  a  maximum  of  5  mature  red  spruce  trees 
for  each  of  44  plots  for  the  years  1900-84.  The 
mean  level  of  residuals  from  forecasts  for  the 
last  20  years  of  the  series  were  generally 
negative,  indicating  a  reduced  ring  width 
relative  to  predicted  ring  width.  These  forecast 
residuals  showed  substantial  spatial  dependence 
that  could  not  be  explained  by  geographical 
factors  alone.  When  both  geographical  and  biotic 
factors,  primarily  measures  of  stand  quality, 
were  taken  into  account,  the  residual  variation 
in  ring  widths  showed  a  weaker  pattern  of  local 
spatial  dependence. 


INTRODUCTION 

The  National  Park  Service  (NPS)  and  the 
Tennessee  Valley  Authority  (TVA)  conducted 
studies  on  experimental  plots  in  the  Great  Smoky 
Mountains,  producing  a  substantial  data  base  of 
information  that  can  be  used  to  examine  annual 
ring  widths  of  red  spruce  (Picea  rubens  Sarg.). 

In  this  report  we  discuss  one  approach  to 
analyzing  ring  widths.  The  objective  of  our 
study  was  to  detect  recent  changes  during  a 
designated  time  series  (1900-84)  that  may  be 
attributable  to  environmental  changes,  such  as 
the  occurrence  of  acid  deposition.  Variations  in 
ring  widths  relative  to  historical  patterns  are 
assessed.  Also  described  is  how  to  determine 
whether  such  patterns  exist  because  of  plot 
characteristics  or  additional  spatial  effects. 

The  main  steps  of  the  study  may  be  summarized 
as  follows: 

1.  Construct  an  average  ring  width  time  series 
for  each  of  the  study  plots  established  by 
the  TVA.  Plots  established  by  the  NPS  were 
not  included. 

2.  Develop  measures  of  recent  increases  or 
decreases  in  ring  width  for  each  plot 
relative  to  forecast  values. 

3.  Relate  the  increases  or  decreases  in  ring 
width  to  geographical  and  biotic  factors. 

4.  Determine  whether  there  is  any  spatial 
pattern  to  the  values  of  excess  or 
deficiency  and  whether  geographical  and 
biotic  factors  are  responsible  for  the 
pattern. 


Step  1. - -Construct  an  average  ring  width  time 
series  for  each  of  the  study  plots  established  by 
the  TVA. 

The  following  are  decisions  made  during  the 
exploratory  stage  of  the  analysis: 

Choice  of  Plots . - -The  analysis  reported  in  this 
project  concerns  48  experimental  plots 
established  by  the  TVA  in  the  Great  Smoky 
Mountains  in  North  Carolina.  The  decision  not  to 
include  the  experimental  plots  established  by  the 
NPS  in  the  same  general  region  was  motivated  by 
time  and  resource  constraints. 

Choice  of  Measurement  Scale. - -Graphs  of  the 
time  series  for  each  of  two  cores  taken  from  five 
trees  usually  produced  very  similar  patterns. 
Since  the  overall  plot  was  the  focus  for  this 
study,  the  two  core  series  were  averaged  for  each 
tree.  The  selection  of  a  designated  time  period 
as  a  series  for  the  entire  plot  was  more 
difficult,  since  individual  trees  may  show 
considerable  variations  from  year  to  year.  The 
overall  mean  was  chosen  as  the  measure  of  average 
ring  width  for  the  plot.  Further  consideration 
of  this  issue  is  presented  in  the  Discussion 
section. 

Selection  of  Trees  and  Time  Frame  for 
Analysis .  --Graphs  for  the  time  series  of  ring 
widths  for  each  of  the  5  trees  per  plot  were  then 
constructed  for  all  48  plots.  Examples  of  four 
of  these  graphs  are  shown  in  figures  1  through  4. 
From  an  examination  of  the  48  graphs,  the 
following  decisions  were  made: 

1 .  Only  red  spruce  would  be  used  in  the  data 
analysis  to  remove  some  heterogeneity  from 
the  time  series  of  ring  widths  averaged 
across  trees  in  a  plot.  Sample  size  was 
not  seriously  reduced  because  214  of  234 
trees  in  the  study  were  red  spruce,  and 
only  4  plots  had  no  red  spruce. 

2.  Only  red  spruce  trees  with  a  pith  date 
earlier  than  1940  would  be  included  in  the 
study.  Therefore  25  red  spruce  trees  were 
eliminated,  and  some  heterogeneity  caused 
by  an  apparent  initial  rapid  increase  in 
ring  width  in  the  early  years  of  growth  was 
alleviated. 

3.  Ring  widths  from  1900  onward  were  analyzed. 
The  heterogeneity  caused  by  the  staggered 
entrance  of  trees  into  the  plot  averages 
and  by  the  apparent  initial  rapid  increase 
in  ring  width  was  minimized. 


J.  Keith  Ord  is  a  professor  in  the  Departments  of  Management  Science  and 
Statistics.  Janice  A.  Derr  is  managing  director  of  the  Statistical 
Consulting  Center,  Pennsylvania  State  University,  University  Park, 
Pennsylvania. 
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Figure  1. - -Graph  of  ring  widths  of  red  spruce  (Picea  rubens  Sarg.)  on  plot 
6  of  48  selected  experimental  plots  established  by  the  Tennessee 
Valley  Authority  in   the  Great  Smoky  Mountains . 
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Figure    2. --Graph  of  ring  widths  of  red  spruce    (Picea  rubens  Sarg.)   on  plot 
18      of      48      selected      experimental      plots      established     by      the 
Tennessee  Valley  Authority  in   the  Great  Smoky  Mountains. 
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Figure   3 . - -Graph  of  ring  widths  of  red  spruce    (Picea  rubsns  Sarg.)   on  plot 
23      of      48      selected      experimental      plots      established     by      the 
Tennessee  Valley  Authority  in   the  Great  Smoky  Mountains . 
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Figure   4. --Graph   of  ring  width   of  red   spruce    (Picea   rubens   Sarg.)    on  plot 

31      of      48      selected      experimental      plots      established      by      the 
Tennessee   Valley  Authority   in    the  Great   Smoky  Mountains . 
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Choice  of  Explanatory  Variables .  - -The  plot 
characteristics  that  were  used  in  the  study 
included  both  geographical  and  biotic  factors  as 
measured  by  survey  teams .  Average  annual 
temperature  and  total  annual  precipitation  for 
the  three  climatic  regions  in  the  study  (North 
Carolina  northern  mountains,  North  Carolina 
southern  mountains,  and  Virginia  southern 
mountains)  appeared  to  be  fairly  similar  in  the 
occurrence  of  peaks  and  dips.  Therefore,  because 
detailed  models  of  precipitation  are  being 
developed  by  others  in  the  project,  climate 
variables  were  not  included  at  this  stage. 

Upon  completion  of  the  exploratory  data 
analysis,  1  time  series  of  ring  widths  was 
constructed  for  each  of  the  44  remaining  plots 


(48  minus  the  4  with  no  red  spruce)  from  the  TVA 
study.  The  time  series  began  with  the  year  1900 
and  ended  with  1984.  The  series  included  only 
those  red  spruce  trees  with  pith  dates  earlier 
than  1940.  Each  time  series  entry  was  an  average 
of  the  width  of  two  cores  from  each  tree ,  taken 
from  a  maximum  of  five  red  spruce  trees .  Table  1 
summarizes  characteristics  of  the  data  for  each 
plot  and  refers  to  geographical  factors,  and 
table  2  refers  to  biotic  factors. 

Step  2. --Develop  measures  of  recent  increases 
or  decreases  in  ring  width  for  each  plot  relative 
to  forecast  values. 

Time  Series  Analysis . - -To  determine  the  recent 
pattern  of  tree  ring  growth  on  each  plot,  an 


Table  l.--PIot  characteristics   and  geographical   variables 


Observation 


Plot 


ELEV 


ASPECT 


LAT 


LONG 


1 

2 

3 

4 

5* 

6 

7 

8 

9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22* 
23 
24 
25 
26 
27 
28 
29 
30* 
31 
32 
33 
34 
35* 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 


1 

5200 

2 

5640 

3 

5200 

4 

6020 

5 

6140 

6 

5700 

7 

5740 

8 

4480 

9 

6000 

10 

6060 

11 

5520 

12 

6200 

13 

6000 

14 

6360 

15 

5880 

16 

6280 

17 

4480 

18 

5470 

19 

5490 

20 

5300 

21 

5340 

22 

5420 

23 

5140 

24 

5100 

25 

5280 

26 

5340 

27 

5300 

28 

4140 

29 

5700 

30 

6015 

31 

4540 

32 

5840 

33 

5040 

34 

5110 

35 

6000 

36 

5780 

37 

4880 

46 

4400 

54 

5120 

55 

5800 

109 

5650 

113 

6240 

125 

5000 

132 

4480 

153 

5060 

158 

4620 

159 

5240 

207 

4880 

21 

36. 

6653 

81, 

5375 

330 

36. 

6597 

81, 

,5472 

210 

36. 

6542 

81. 

,5333 

35 

36. 

1055 

82 

,1333 

213 

36 

1000 

82 

,1222 

90 

35 

3278 

82 

,9612 

16 

35 

,3500 

82 

,9612 

205 

36 

.3417 

81 

,6500 

220 

35 

.7375 

82 

,3195 

225 

35 

,7250 

82 

,2917 

220 

35 

.7292 

82 

,2792 

245 

35 

,7445 

82 

,3250 

255 

35. 

8305 

82 

,2555 

137 

35 

.7750 

82 

,2583 

303 

35 

8167 

82 

,2555 

300 

35 

,7833 

82 

,2583 

340 

36 

,3500 

81 

.6417 

190 

36. 

6375 

81 

.6055 

345 

36. 

,6403 

81 

.6055 

55 

36. 

6638 

81 

5403 

55 

35. 

7208 

82 

.2792 

142 

35. 

7458 

82 

.2722 

125 

35. 

7417 

82 

.2750 

190 

35. 

7458 

82 

2667 

305 

35. 

,7333 

82 

.3250 

123 

35 

,7250 

82 

,3125 

301 

35 

,7222 

82 

,3083 

23 

35 

,7250 

82 

,2722 

30 

36 

,1042 

81 

.8125 

345 

36 

,0917 

82 

.1500 

235 

36 

.3388 

81 

.6542 

0 

35 

.4638 

83 

.1375 

90 

36 

,1112 

81 

.7958 

290 

36 

.1083 

81 

.8250 

305 

35 

.8250 

82 

.2555 

120 

35 

.8458 

82 

.2417 

193 

35 

.8388 

82 

.2388 

50 

36 

.3458 

81 

.6500 

35 

35 

.4750 

83 

.1167 

268 

35 

.4750 

83 

.0958 

220 

35 

.7888 

82 

.2667 

265 

35 

.7292 

82 

.2917 

350 

36 

.1083 

82 

.1333 

215 

36 

.1292 

82 

.2917 

305 

35 

.2945 

82 

.9375 

20 

35 

.3625 

82 

.8833 

60 

35 

.3388 

82 

.9612 

340 

35 

,3583 

82 

,8542 
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Plots  eliminated  from  study  because  no  red  spruce  was  present. 


Table  2. --Plot   characteristics   and   biotic   variables 


Observation 

SBALIVE 

SBADEAD 

SB  AX 

LIVETREE 

DEADTREE 

TREEX 

1 

56.7 

4.3 

0.92951 

6672 

185 

0.97302 

2 

38.6 

5.9 

0.86742 

10329 

247 

0.97665 

3 

39.6 

7.3 

0.84435 

7771 

1186 

0.86759 

4 

11.6 

11.1 

0.51101 

642 

1482 

0.30226 

5* 

41.5 

6.0 

0.87368 

3064 

544 

0.84922 

6 

30.2 

3.3 

0.90149 

4695 

704 

0.86961 

7 

31.2 

3.7 

0.89398 

1507 

408 

0.78695 

8 

43.1 

0.5 

0.98853 

1742 

148 

0.92169 

9 

48.2 

3.5 

0.93230 

9798 

111 

0.98880 

10 

37.5 

6.3 

0.85616 

4547 

111 

0.97617 

11 

30.1 

2.8 

0.91489 

1532 

988 

0.60794 

12 

41.7 

26.6 

0.61054 

10687 

185 

0.98298 

13 

14.7 

2.6 

0.84971 

1705 

309 

0.84657 

14 

26.6 

7.2 

0.78698 

1915 

741 

0.72101 

15 

25.1 

12.5 

0.66755 

4176 

1161 

0.78246 

16 

12.0 

12.6 

0.48780 

7141 

543 

0.92933 

17 

39.0 

2.4 

0.94203 

2508 

222 

0.91868 

18 

55.7 

3.3 

0.94407 

1174 

383 

0.75401 

19 

34.2 

1.5 

0.95798 

1507 

74 

0.95319 

20 

45.1 

4.9 

0.90200 

2002 

99 

0.95288 

21 

45.9 

1.3 

0.97246 

1149 

173 

0.86914 

22* 

52.8 

2.6 

0.95307 

8698 

124 

0.98594 

23 

37.7 

0.0 

1.00000 

531 

0 

1.00000 

24 

60.7 

1.1 

0.98220 

3286 

222 

0.93672 

25 

43.5 

19.6 

0.68938 

3632 

136 

0.96391 

26 

69.0 

0.0 

1.00000 

7401 

0 

1.00000 

27 

52.5 

8.8 

0.85644 

3410 

432 

0.88756 

28 

63.2 

10.1 

0.86221 

840 

210 

0.80000 

29 

18.5 

15.1 

0.55060 

4707 

334 

0.93374 

30* 

45.0 

10.0 

0.81818 

9044 

2557 

0.77959 

31 

35.4 

2.7 

0.92913 

1680 

297 

0.84977 

32 

10.5 

32.4 

0.25060 

1075 

1631 

0.39727 

33 

18.9 

2.9 

0.86697 

1025 

86 

0.92259 

34 

36.3 

4.9 

0.88107 

1248 

556 

0.69180 

35* 

15.5 

0.2 

0.98726 

1520 

25 

0.98382 

36 

38.7 

13.7 

0.73855 

4213 

272 

0.93935 

37 

59.6 

0.9 

0.98512 

914 

12 

0.98704 

38 

34.1 

0.5 

0.98555 

1124 

148 

0.88365 

39 

5.9 

0.8 

0.88060 

37 

80 

0.31624 

40 

39.5 

2.4 

0.94272 

1606 

219 

0.88000 

41 

28.2 

4.9 

0.85196 

2718 

292 

0.90299 

42 

15.3 

0.0 

1.00000 

3497 

0 

1.00000 

43 

38.4 

4.1 

0.90353 

8686 

1819 

0.82684 

44 

46.0 

0.4 

0.99138 

1594 

70 

0.95793 

45 

46.7 

0.9 

0.98109 

1557 

173 

0.90000 

46 

21.6 

1.9 

0.91915 

1890 

400 

0.82533 

47 

10.5 

7.8 

0.57377 

1408 

2088 

0.49275 

48 

28.2 

19.1 

0.59619 

803 

738 

0.52109 

Plots  eliminated  from  study  because  no  red  spruce  was  present. 


autoregressive-integrated-moving  average  (ARIMA) 
time  series  model  for  each  of  the  plots  was  first 
developed.  Kendall  and  others  (1983)  and 
Vandaele  (1983)  provide  details  of  ARIMA  models 
and  the  underlying  assumptions  . 

Time  Series  Modeling . - -Vafious  years  of  the 
series  indicated  marked  trends  of  tree  ring 
growth.  To  accommodate  trends,  the  series  were 
differenced  where  necessary.  Other  approaches  to 
this  problem  are  covered  in  the  Discussion 
section.  In  this  data  set,  it  was  never 
necessary  to  difference  more  than  once.  A 
summary  of  the  fitted  models  is  presented  in 
table  3.   From  the  table,  it  can  be  seen  that 


many  of  the  series  were  described  satisfactorily 
by  an  autoregressive  scheme  of  order  2, 
occasionally  with  higher  order  moving  average 
(MA)  terms. 

The  MA  terms  improved  the  fit  as  measured  by 
the  diagnostics  but  did  not  materially  affect  the 
forecasts.  The  autoregression  [AR^2>] 
coefficients  were  usually  both  positive  with  <j>i  + 
4>2  in  the  range  0.6  to  0.9,  indicating  a  carry- 
over from  one  growing  season  %p  the  next,  as 
would  be  expected.  When  <j>i  +  <j>2  exceeded  0.9, 
nonstationarity  in  the  series  was  evident,  and 
differencing  was  performed.  A  low-order  MA 
scheme  usually  gave  an  adequate  description  of 
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Table  3 .- -Summary   of   autoregressive- integrated-moving   average   models   for 
each  plot 


Plot  No. 

Difference 

Lag  Structure"4" 

Comments 

AR 

MA 

1 

No 

2 

0 

2 

No 

2 

0 

3 

No 

2 

0 

4 

No 

2 

0 

6 

No 

2 

0 

7 

No 

2 

(7?) 

8 

Yes 

0 

(1,4,10) 

Short  series 

9 

No 

2 

0 

10 

No 

2 

0 

11 

No 

3 

0 

Short  series 

12 

No 

2 

0 

Short  series 

13 

Yes 

0 

1 

Short  series 

14 

Yes 

2 

0 

Short  series 

15 

Yes 

(3) 

0 

Short  series 

16 

Yes 

0 

(4) 

Short  series 

17 

Yes 

0 

2 

18 

No 

2 

(8) 

19 

No 

2 

0 

20 

Yes 

0 

2 

21 

Yes 

1 

0 

23 

Yes 

0 

1 

24 

Yes 

0 

3 

25 

No 

2 

0 

26 

No 

2 

0 

27 

No 

2 

0 

28 

No 

2 

0 

29 

Yes 

0 

(2,9) 

31 

Yes 

0 

(1,2,6) 

32 

No 

2 

(6) 

33 

Yes 

0 

(1,2,6) 

34 

No 

2 

0 

36 

No 

2 

0 

37 

No 

2 

(10?) 

Structural 
change  in 

46 

No 

1 

0 

series? 

54 

No 

2 

0 

55 

No 

2 

0 

109 

Yes 

0 

2 

First  seven 
terms  deleted 

113 

Yes 

0 

0 

125 

Yes 

0 

2 

Short  series 

132 

Yes 

0 

1 

153 

Yes 

0 

(1,4,5,6) 

158 

Yes 

0 

(2,5) 

159 

Yes 

0 

(1,4) 

Short  series 

207 

No 

2 

(3,7?) 

k  indicates  lags  1,  2 k 

(k)  indicates  lag  k  only 
(j,k)  indicates  lags  j,k  only 
(k?)  indicates  lag  k  a  possibility 
AR  =  autoregression  coefficient 
MA  =  moving  average 
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the  differenced  series.  It  is  well  known  that  AR 
schemes  with  a  root  of  the  auxiliary  equation 
near  unity  can  often  be  well  approximated  by  a 
low-order  MA  scheme  with  a  single  difference. 
Therefore  the  models  are  not  very  different  in 
practice  despite  their  distinct  theoretical 
properties . 

The  only  series  that  caused  major  problems  was 
that  for  plot  109,  where  major  increases  in  the 
first  7  years  were  followed  by  steady  declines. 
After  the  data  for  the  first  7  years  were 
deleted,  a  satisfactory  model  was  fitted.  It  can 
be  assumed  that  the  use  of  ring  width  rather  than 
incremental  basal  area  was  the  cause  of  these 
nonstationarity  problems. 

Measures  of  Recent  Relative  Change  in  Ring 
Width. - -To  assess  recent  relative  increases  or 
decreases  in  ring  width,  each  series  over  the 
periods  was  forecasted: 

1.  1965-1984,   using   1964   as   the   forecast 
origin. 

2.  1975-1984,   using   1974   as   the   forecast 
origin. 

The  residuals,  the  difference  -between  observed 
and  predicted  values,  were  then  computed  for 
each  year  in  the  period.  The  means  and  standard 
deviations  of  these  residuals  were  computed  for 
each  plot  (table  4).  It  should  be  noted  that  the 
models  were  fitted  to  the  entire  series,  1900-84, 
and  forecasts  were  then  generated  from  the 
forecast  origin.   A  pure  forecasting  method  would 


have  involved  fitting  to  1964  (or  1974)  and  then 
forecasting.  However,  the  risk  of  structural 
changes  in  the  series  was  such  that  the  pure 
forecasts  might  misrepresent  recent  trends. 
Although  our  approach  biases  the  residuals 
somewhat  towards  zero,  the  method  seemed  to 
provide  a  clearer  picture  of  recent  developments. 
Because  changes  in  ring  width  might  be 
considered  in  either  absolute  or  percentage 
terms,  also  considered  were  the  indicators: 


proportional  change  = 


average  of  residuals 
average  ring  width 


over  the  two  forecast  periods.   These  values  are 
also  listed  in  table  4  as  PCT20  and  PCT10 . 

Assessment  of  Mean  Change .  - -One  should  note 
whether  the  residual  ring  widths  are  below  the 
expected  value  of  zero  for  the  periods 
considered.  The  results  of  one-tailed  t-tests  on 
the  data  in  table  4  were  as  follows: 


Variable 


Mean 


t-value 


Adjusted 
t-value 


RES20 

-17.47 

-3.28 

-1.89 

PCT20 

-0.108 

-3.46 

-2.34 

RES  10 

4.43 

0.93 

0.50 

PCT10 

0.022 

0.79 

0.47 

Table  4 .-- Summary  statistics  from   time   series   analysis 


Summarv  Statistics^ 

Observation 

OVMEAN 

RES20 

R20SD 

PCT20 

RES10 

R10SD 

PCT10 

1 

189.5 

18.3 

34.4 

0.9657 

-14.3 

43.3 

-0.07546 

2 

158.5 

-34.5 

25.3 

-0.21767 

-4.8 

21.5 

-0.03028 

3 

187.3 

-48.1 

24.8 

-0.25681 

-1.0 

9.9 

-0.00534 

4 

5* 

6 

147.1 

-50.5 

20.6 

-0.34330 

-15.6 

22.1 

-0.10605 

198.6 

-3.7 

30.4 

-0.01863 

5.2 

33.7 

0.02613 

7 

147.9 

-37.1 

23.1 

-0.25085 

-35.1 

26.9 

-0.23732 

8 

375.3 

44.8 

63.5 

0.11937 

93.7 

63.9 

0.24967 

9 

140.0 

-42.1 

30.8 

-0.30071 

-22.1 

16.0 

-0.15786 

10 

138.7 

-31.5 

23.9 

-0.22711 

-28.8 

29.9 

-0.20764 

11 

321.9 

-57.7 

52.1 

-0.17925 

-30.2 

43.0 

-0.09382 

12 

163.6 

-16.3 

22.4 

-0.09963 

-5.4 

26.9 

-0.03301 

13 

152.0 

17.4 

34.7 

0.11447 

62.6 

35.6 

0.41184 

14 

259.1 

25.8 

31.4 

0.09958 

41.1 

31.3 

0.15863 

15 

132.4 

-18.5 

18.6 

-0.13973 

26.6 

21.8 

0.20091 

16 

134.9 

84.4 

75.5 

0.62565 

85.6 

48.1 

0.63454 

17 

189.9 

3.4 

30.6 

0.01790 

15.4 

13.7 

0.08110 

18 

182.3 

-39.0 

22.5 

-0.21393 

-17.6 

29.1 

-0.09654 

19 

81.6 

-32.7 

8.3 

-0.40074 

6.6 

7.6 

0.08088 

20 

232.5 

-58.1 

37.1 

-0.24989 

-9.2 

37.3 

-0.03957 

21 
22* 

239.0 

-113.1 

48.9 

-0.47322 

-39.7 

43.5 

-0.16611 

23 

182.4 

-25.1 

36.9 

-0.13761 

24.3 

48.5 

0.13322 

24 

192.5 

-45.8 

30.4 

-0.23792 

36.6 

32.0 

0.19013 

25 

126.5 

-26.4 

21.0 

-0.20870 

-16.1 

24.8 

-0.12727 

26 

172.6 

-20.7 

31.2 

-0.11993 

-23.5 

40.5 

-0.13615 

27 

131.2 

-31.2 

28.2 

-0.23780 

-24.0 

32.3 

-0.18293 

28 

125.5 

-19.3 

16.2 

-0.15378 

-10.3 

19.5 

-0.08270 

29 
30* 

82.7 

8.4 

22.2 

0.10157 

27.7 

23.2 

0.33495 

31 

203.9 

-55.5 

20.4 

-0.27219 

-5.8 

22.1 

-0.02845 

32 

165.2 

-15.9 

24.2 

-0.09625 

31.6 

26.0 

0.19128 

33 

184.3 

-3.4 

37.7 

-0.01845 

26.6 

48.5 

0.14433 

34 

173.8 

21.3 

30.8 

0.12255 

26.8 

38.9 

0.15420 

35' 
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Table  4 .-- Summary  statistics  from   time   series   analysis- -Continued 


Summa 

rv  Statistic 

s± 

Observation 

OVMEAN 

RES20 

R20SD 

PCT20 

RES10 

R10SD 

PCT10 

36 

103.1 

-19.4 

16.5 

-0.18817 

-10.1 

23.5 

-0.0979 

37 

112.1 

-12.4 

15.8 

-0.11062 

-5.8 

21.7 

-0.05174 

38 

189.5 

-14.5 

41.8 

-0.07652 

-24.5 

51.2 

-0.12929 

39 

136.5 

-27.4 

20.9 

-0.20073 

-21.9 

25.1 

-0.16044 

40 

236.2 

-12.5 

31.0 

-0.05292 

-20.8 

33.1 

-0.08806 

41 

288.5 

41.0 

31.6 

0.14211 

33.4 

21.3 

0.11577 

42 

187.7 

-22.1 

32.3 

-0.11774 

-3.5 

35.6 

-0.01865 

43 

185.9 

24.6 

36.8 

0.13233 

50.7 

39.6 

0.27273 

44 

189.2 

-4.9 

33.4 

-0.02590 

23.6 

39.9 

0.12474 

45 

139.3 

-88.8 

34.8 

-0.63747 

-35.7 

30.4 

-0.25628 

46 

176.9 

1.3 

23.5 

0.00735 

-7.8 

32.2 

-0.04409 

47 

273.3 

-45.2 

37.5 

-0.16539 

-7.2 

37.6 

-0.02634 

48 

127.9 

14.2 

13.7 

0.11102 

17.5 

14.8 

0.13683 

OVMEAN 


overall  mean  of  series. 


RES20  =  residuals  from  forecasts  for  last  20  years. 

R20SD  =  standard  deviation  of  RES20  values. 

PCT20  =  RES20/OVMEAN. 

RES10,  R10SD,  PCT10  are  defined  similarly. 

Plots  eliminated  from  study  because  no  red  spruce  was  present. 


The  adjusted  t-values  were  computed  following  the 
approach  described  by  Cliff  and  Ord  (1981) 
modified  to  the  one -sample  case.  The  adjustment 
takes  account  of  the  positive  spatial  dependence 
among  the  data  and  may  be  written  as: 

tadj  =  t  (1  -  I), 

where  I  is  defined  in  equation  (1)  under  Step  4. 

Evidently  the  20-year  residual  ring  widths 
have  a  mean  that  is  significantly  less  than  zero, 
while  the  null  hypothesis  of  a  zero  mean  is 
accepted  for  the  10-year  values.  Therefore  a 
drop  is  indicated  in  average  ring  width  in  the 
1960 's  that  has  subsequently  stabilized  at  that 
lower  level . 

Step  3. --Relate  recent  increases  and  decreases 
in  ring  width  to  geographical  and  biotic  factors. 

In  this  step,  the  changes  in  ring  width  were 
related  to  the  various  plot  characteristics  to 
determine  if  there  were  any  explanation  for  the 
changes . 

Regression  Analysis .  - -Each  of  the  four  residual 
ring  width  measures  was  modeled  using  stepwise 
regression  with  the  following  variables: 

geographical:  latitude  (LAT) ,  longitude 
(LONG),  (latitude)2  -  LAT2 ,  (longitude)2  = 
LONG 2,  latitude  *  longitude  =  LATLONG, 
elevation  (ELEV) ,  and  aspect  (coded  as  sine 
and  cosine,  SASP  and  CASP) . 
biotic:  number  of  live  trees  (LIVETREE) , 
number  of  dead  trees  (DEADTREE) ,  stand  basal 
area  of  live  trees  (SBALIVE) ,  stand  basal 
area  of  dead  trees  (SBADEAD) ,  and  two  derived 
indices : 


TREEX  = 
SBAX  = 


LIVETREE 


(LIVETREE  +  DEADTREE) 

SBALIVE 

(SBALIVE  +  SBADEAD) 


where  TREEX  is  proportion  of  live  trees,  and  SBAX 
is  the  proportion  of  live  basal  area.  The  values 
of  these  variables  are  listed  in  table  1. 

The  quadratic  factors  of  latitude  and  longitude 
were  included  to  allow  a  low- order  trend  surface 
analysis  (Cliff  and  others  1975).  However, 
initial  runs  using  only  the  geographical 
variables  showed  virtually  no  correlation  between 
any  of  these  variables  and  the  residual  ring 
width  measures ;  therefore  they  have  not  been 
reported  separately.  A  total  of  eight  analyses 
are  reported  in  tables  5  through  8.  For  each  of 
the  residual  ring  width  measures,  the  analysis 
was  performed  using  both  unweighted  and  weighted 
least  squares  (LS)  .  The  weights  used  were  the 
standard  deviations  given  in  table  4.  In  all 
cases,  the  significance  level  for  a  variable  to 
leave  or  stay  was  set  at  0.25.  The  residuals 
from  the  regression  analyses  are  given  in  tables 
9  and  10. 

Interpretation  of  Regression  Results.  -- 
Comparisons  within  and  across  tables  5  through  8 
show  the  following: 

1.  The  value  of  R2  is  in  all  cases  somewhat 
higher  for  weighted  LS  than  for 
unweighted.  A  high  standard  deviation  in 
the  time  series  residuals  shows  an  erratic 
ring  width  pattern.  Therefore  the 
weighting  is  useful  because  greater 
emphasis  is  given  to  the  plots  with  more 
stable  ring  width  development.  Otherwise 
the  same  variables  were  selected  by  the 
stepwise  procedures  for  both  estimation 
procedures,  and  the  two  sets  of  estimates 
were  broadly  consistent  for  each  of  the 
four  dependent  variables. 

The  proportional  change  indicators  yield 
models  with  a  higher  degree  of  explanatory 
power  than  those  based  on  absolute  changes . 
Since  average  ring  widths  vary  considerably 
between  sites,  use  of  the  proportional 
change  indicator  seems  preferable . 


2. 
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Table  5 . - -Regression   analysis  for   Che   dependent   variable  20-year  mean 
residual   ring  width    (RES20) 


Unweighted  least  squares 


Sot 


df 


Sum  of 
squares 


Mean 
square 


F  value 


Prob>F 


Model 

4 

10813.450 

2703.362 

Error 

39 

42758.529 

1096.373 

C  total 

43 

53571.979 

Root  MSE 

33.111517 

R- square 

Dep  mean 

-17.465909 

Adj  R-sq 

C.V. 

-189.578 

2.466    0.0608 


0.2018 
0.1200 


Variable 


df 


Parameter 
estimate 


Standard 
error 


T  for  HO: 
Parameter=0 


Prob  >  |T| 


Intercept 

1 

70.597735 

70.426505 

1.002 

0.3223 

ELEV 

1 

-0.013709 

0.010648 

-1.287 

0.2055 

SBAX 

1 

-69.517093 

40.250901 

-1.727 

0.0921 

TREEX 

1 

91.734416 

36.388501 

2.521 

0.0159 

SBALIVE 

1 

-0.916871 

0.442469 

-2.072 

0 . 0449 

Weighted  least  squares 


Sum  of 

Mean 

Source 

df 

squares 

square 

F  value 

Prob>F 

Model 

4 

800896 

200224 

4.859 

0.0028 

Error 

39 

1607051 

41206.433 

C  total 

43 

2407947 

Root  MSE 

202.994 

R- square 

0.3326 

Dep  mean 

-13.665677 

Adj  R-SQ 

0.2642 

C.V. 

-1485.43 

Parameter 

Standard 

T  for  HO: 

Variable 

df 

estimate 

error 

Parameter=0 

Prob  >  |T| 

Intercept 

1 

99.:87126 

81.765087 

1.213 

0.2324 

ELEV 

1 

-0.017920 

0.011910 

-1.505 

0.1405 

SBAX 

1 

-105.865 

45.857356 

-2.309 

0.0264 

TREEX 

1 

135.208 

39.969084 

3.408 

0.0015 

SBALIVE 

1 

-1.201380 

0.502896 

-2.389 

0.0218 

The  regression  analyses  for  geographical 
variables  provided  only  very  little 
explanatory  power.  When  the  biotic 
variables  were  also  included,  elevation 
became  important,  and  the  residual  ring 
width  became  more  negative  as  elevation 
increased.  This  suggests  that  the  higher 
elevation  plots  did  worse  than  average  over 
the  10-  and  20-year  periods  considered. 
The  only  other  geographical  variables  that 
appeared  in  any  models  were  L0NG2  and  SIN 
(aspect) .  The  coefficient  on  L0NG2 
suggests  a  downward  trend  in  the  10-year 
change  variables  from  east  to  west.  Since 
the  plot  locations  extended  approximately 
northeast  to  southwest,  this  may  reflect 
the  influence  of  climatic  factors.  The 
coefficient  for  SIN  (aspect)  indicates  that 
the  proportional  change  variable  for  the 
10-year  period  is  higher  in  plots  with  a 
southerly  aspect.  Again,  this  may  reflect 
climatic   effects. 


4.  The  most  important  variable  in  almost  all 
cases  was  the  tree  index,  TREEX,  which  is 
probably  an  indicator  of  stand  health,  and 
strong,  positive  correlation  is  to  be 
expected.  The  other  major  biotic  variable 
was  SBALIVE,  but  this  appears  with  a 
negative  sign  in  the  regression.  SBADEAD 
and  the  stand  basal  area  index  (SBA)  also 
appear  on  occasion,  again  with  negative 
signs  in  all  cases.  The  interpretation  of 
these  effects  is  unclear,  but  these 
variables  may  relate  to  other  biotic 
factors  such  as  the  age  of  the  stand  and 
the  degree  of  competition. 

Overall,  the  weighted  regressions  on  the 
proportional  change  indicators  appear  to  give 
a  reasonable  explanation  of  the  variations  in 
residual  ring  width. 

Step  4. --Determine  whether  there  is  any  spatial 
pattern  to  the  values  of  excess  or  deficiency  and 
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Table  6 . - -Regression   analysis  for  10-year  mean   residual   ring  width    (RES10) 


Unweighted  least  squares 


Source 


df 


Sum  of 
squares 


Mean 
square 


F  value 


Prob>F 


Model 

2 

6777.608 

3388.804 

3 

.939 

Error 

41 

35269.460 

860.231 

C  total 

43 

42047.067 

Root  MSE 

29.329690 

R- square 

0 

.1612 

Dep  mean 

4.427273 

Adj  R-sq 

0 

.1203 

C.V. 

662.4776 

0.0272 


Variable 


df 


Parameter 
estimate 


Standard 
error 


T  for  HO: 
Parameter=0 


Prob  >  | T | 


Intercept 

1 

814.995 

404.974 

2.012 

0.0508 

SBALIVE 

1 

-0.763081 

0.303229 

-2.517 

0.0159 

LONG  2 

1 

-0.115851 

0.059378 

-1.951 

0.0579 

Weighted  least  squares 


Sum  of 

Mean 

Source 

df 

squares 

square 

F  value 

Prob>F 

Model 

2 

301627 

150813 

4.325 

0.0198 

Error 

41 

1429655 

34869.627 

C  total 

43 

1731281 

Root  MSE 

186.734 

R- square 

0.1742 

Dep  mean 

7.603876 

Adj  R-sq 

0.1339 

C.V. 

2455.775 

Parameter 

Standard 

T 

for  HO: 

Variable 

df 

estimate 

error 

parameter=0 

Prob  >  |T| 

Intercept 

1 

1097.709 

474.049 

2.316 

0.0257 

SBALIVE 

1 

-0.852323 

0.342500 

-2.489 

0.0170 

LONG  2 

1 

-0.156754 

0.069541 

-2.254 

0.0296 

whether  this  can  be  accounted  for  by 
geographical  and  biotic  factors . 

In  this  section,  the  spatial  methods  used  are 
first  described.  Then  the  spatial  analysis  for 
the  residual  ring  widths  and  for  their  residuals 
from  the  regression  equations  developed  in  step  3 
is  presented.  The  objective  of  the  spatial 
analysis  is  to  discover  if  there  is  any  spatial 
pattern  in  the  recent  changes  in  ring  width,  both 
among  the  initial  values  and  the  residuals,  from 
the  regression  equations. 

Spatial  Methods .  - -The  first  step  in  any  spatial 
analysis  is  to  determine  whether  or  not  there  is 
any  evidence  of  spatial  pattern  among  the  data, 
given  the  plot  locations.   If  the  n  plots  have 


observed  values  x^  (i  =  1, 


,  n) ,  we  set  Zj 


and  use  the   spatial  autocorrelation 


statistic : 


Under  the  null  hypothesis  (Hq)  of  no  spatial 
autocorrelation  (or  independence) ,  it  may  be 
shown  that: 


E(I)  -  -l/(n  -1). 

Cliff  and  Ord  (1981)  show  the  distribution  of  I 
under  Hq  to  be  approximately  normal,  provided 
that  n  is  not  too  small .  For  the  configurations 
of  weights  used  and  the  number  of  plots  available 
(n  =  44)  ,  the  normal  approximation  is 
satisfactory.  It  should  be  noted  that  I  is  not 
defined  quite  like  a  regular  correlation 
coefficient;  in  particular,  the  values  tend  to  be 
closer  to  the  origin  than  one  would  expect.  For 
this  reason,  the  magnitudes  of  the  standard 
deviates 


n  II  WU  zi  zj 


I  = 


i  i 


So  I  < 


i  J 


(1) 


where  S0  =  ^  I  wij  and  the  (wi,  are  a  set  of 


negative  weights  to  be  specified,  with  w±i 


[I  -  E(I)]/4vTTy 

are  often  more  useful  than  the  values  of  I 
themselves.  From  Cliff  and  Ord  (1981),  the 
variance  of  I  under  Hq  is : 

V(I)  =  E(I2)  -  [E(I)]2, 
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Table  7 . - -Regression   analysis   for  20-year  mean   residual   ring  width 
divided  by  overall  mean  ring  width    (PCT20) 


Unweighted 

least  squar 

es 

Sum  of 

Mean 

Source 

df 

squares 

square 

f  value 

Prob>F 

Model 

5 

0.466992 

0.093398 

2.592 

0.0411 

Error 

38 

1.369366 

0.036036 

C  total 

43 

1.836359 

Root  MSE 

0.189831 

R- square 

0.2543 

Dep  mean 

-0.107706 

Adj  R-sq 

0.1562 

C.  V. 

-176.249 

Parameter 

Standard 

T  for  HO: 

Variable 

df 

estimate 

error 

parameter=0 

Prob  >  |T| 

Intercept 

1 

0.995842 

0.583306 

1.707 

0.0959 

ELEV 

1  • 

-.0000856071 

.00006109282 

-1.401 

0.1692 

SBAX 

1 

-1.074991 

0.566563 

-1.897 

0.0654 

SBADEAD 

1 

-0.013608 

0.011104 

-1.225 

0.2279 

TREEX 

1 

0.579832 

0.215236 

2.694 

0.0105 

SBALIVE 

1 

-0.00389848 

0 

.002933337 

-1.329 

0.1918 

Weighted  least  squares 


Sum  of 

Mean 

Source 

df 

squares 

square 

F  value 

Prob>F 

Model 

5 

35.672534 

7.134507 

6.069 

0.0003 

Error 

38 

44.671739 

1.175572 

C  total 

43 

80.344273 

Root  MSE 

1.084238 

R- square 

0 . 4440 

Dep  mean 

-0.073059 

Adj  R-sq 

0.3708 

C.  V. 

-1484.05 

Parameter 

Standard 

T  for  HO: 

Variable 

df 

estimate 

error 

Parameter=0 

Prob  >  |T| 

Intercept 

1 

1.311741 

0.583284 

2.249 

0.0304 

ELEV 

1 

-.0000828751 

.00006361611 

-1.303 

0.2005 

SBAX 

1 

-1.582828 

0.531610 

-2.977 

0.0050 

SBADEAD 

1 

-0.022603 

0.011481 

-1.969 

0.0563 

TREEX 

1 

0.800385 

0.214337 

3.734 

0.0006 

SBALIVE 

1 

-0.00426758 

0.003139655 

-1.359 

0.1821 

where 

E(I2)  =  (n[(n2  -  3n  +  3)SX  -  n  S2  +  3Sq] 
b2[(n2  -  n)Sx  -  2n  S2  +  6S2,]  }/ 
(n  -  1)  (n  -  2)  (n  -  3)Sj}, 

where    n   -  number  of  plots, 

S0  -  I  I   wij  • 

i  J 

Si  -  I   X  (wfj  +  Wij  w  )  , 

i  J 
S2   =  I  (Wi.  +  W.i)2'  and 


i  J 

I 

i 

I 
j 


I   Wy  and  w_d  =  I   Wjl. 


(2) 


Choice  of  Weights .  - -Given  the  irregular  array 
of  plot  locations,  the  choice  of  weights  for  use 
in  (1)  is  somewhat  arbitrary.  However,  when  the 
variables  X  are  normally  distributed,  it  is 
known  (Cliff  and  Ord  1981)  that  I  is  the  locally 
most  powerful  test  for  alternatives  of  the  form: 


HX:   Var(X)  -  a2 [I  -  SW]"1 


(3) 


where  &  and  ol  are  parameters  and  W  is  symmetric. 
The  null  hypothesis  then  becomes  Hq:  &  -  0,  or 
Var(X)  =  a       I.    Since  we  are  interested  in 


detecting   local 
suggests  that  w^^ 


close,   but 


'ij 


spatial  similarities,  this 
>  0  when  plots  i  and  j  are 
0   when   they   are   distant. 
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Table  8 . - -Regression   analysis  for  10-year  mean   residual   ring  width 
divided  by  overall  mean  ring  width    (PCT10) 


Unweighted 

least  squares 

Sum  of 

Mean 

Source 

df 

squares 

square 

F  value 

Prob>F 

Model 

7 

0.678485 

0.096926 

4.467 

0.0012 

Error 

36 

0.781208 

0.021700 

C  total 

43 

1.459693 

Root  MSE 

0.147310 

R- square 

0.4648 

Dep  mean 

0.021891 

Adj  R-sq 

0.3607 

C.  V. 

672.9362 

Parameter 

Standard 

T  for  HO: 

Variable 

df 

estimate 

error 

Parameter=0 

Prob  >  |T| 

Intercept 

1 

3.999241 

2.193178 

1.823 

0.0765 

ELEV 

1  - 

.0000812191 

.00004886896 

-1.662 

0.1052 

SASP 

1 

-0.053525 

0.035049 

-1.527 

0.1355 

LONG  2 

1  - 

0.000375621 

0.0003257647 

-1.153 

0.2565 

SBALIVE 

1 

-0.00354181 

0. 

002301411 

-1.539 

0.1326 

SB  AX 

1 

-1.411441 

0.455795 

-3.097 

0.0038 

SBADEAD 

1 

-0.020464 

0. 

009001847 

-2.273 

0.0291 

TREEX 

1 

0.520329 

0.181001 

2.875 

0.0067 

Weighted  least  squares 


Sum  of 

Mean 

Source 

df 

squares 

square 

F  value 

Pr'ob>F 

Model 

7 

27.779442 

3.968492 

5.605 

0.0002 

Error 

36 

25.489160 

0.708032 

C  total 

43 

53.268603 

Root  MSE 

0.841447 

R- square 

0.5215 

Dep  mean 

0.035922 

Adj  R-sq 

0.4285 

C.V. 

2342.437 

Parameter 

Standard 

T  for  HO: 

Variable 

df 

estimate 

error 

Parameter=0 

Prob  >  |T| 

Intercept 

1 

4.831085 

2.305746 

2.095 

0.0432 

ELEV 

1 

-0000758574 

.00005072042 

-1.496 

0.1435 

SASP 

1 

-0.069365 

0.033810 

-2.052 

0.0475 

L0NG2 

1 

-0.000484338 

0.0003417928 

-1.417 

0.1651 

SBALIVE 

1 

-0.00324435 

0.002267685 

-1.431 

0.1611 

SBAX 

1 

-1.564668 

0.436777 

-3.582 

0.0010 

SBADEAD 

1 

-0.024429 

0.009060183 

-2.696 

0.0106 

TREEX 

1 

0.554004 

0.181281 

3.056 

0.0042 

Furthermore ,  given  the  rapid  changes  in  local 
topography  that  are  possible  in  the  study  area, 
it  is  reasonable  to  set  Wji  =0  when  there  are 
several  plots  between  i  and  j .  Given  that  (1)  is 
being  used  primarily  as  an  exploratory  device, 
these  guidelines  may  be  incorporated  into  the  set 
of  weights  by  use  of  nearest  neighbor  linkages 
(Cliff  and  others  1975).  The  exact  specification 
of  weights  is  not  critical.  Thus  two  sets  of 
weights  are  considered: 


2. 


=  1,  if  plots  i  and  j  are  first  or 


'ij 


—  1,  if  plots  i  and  j  are  nearest 

neighbors,  otherwise 
=  0. 


second  nearest  neighbors,  otherwise 
=  0. 
These  weights  are  not  symmetric,  but  this  does 
not  cause  any  problems.  In  each  case,  a  distance 
criterion  was  used  to  eliminate  linkages  across 
very  long  distances.  The  sets  of  neighbors  are 
summarized  in  (B)  of  the  Appendix.  These  weights 
were  used  in  all  subsequent  analyses.  A  program 
listing  for  a  simple  FORTRAN  program  to  compute  I 
and  the  corresponding  standard  deviate  is  listed 
in  the  Appendix  (A)  .  For  the  first  nearest 
neighbor,  Sq  =  42 ,  S^  =  68,  and  S2  =  192.   For 
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Table  9 . - -Residuals   from  ordinary   least   squares   regression  for  RES20 , 
RES 10,    PCT20,    PCT10,    respectively 


Observation 

RRES20 

RRES10 

PCRES20 

PCRES10 

1 

46.334 

-15.807 

0.26047 

-0.04145 

2 

-21.678 

-19.937 

-0.13374 

-0.12678 

3 

-31.993 

-15.635 

-0.14916 

-0.08203 

4 

5* 

6 

-20.136 

-40.225 

-0.25344 

-0.21953 

14.431 

10.602 

0.10100 

0.12306 

7 

-10.443 

-28.935 

-0.07859 

-0.13183 

8 

59.306 

83.942 

0.21011 

0.17438 

9 

-12.145 

-15.249 

-0.11850 

-0.11506 

10 

-14.668 

-30.645 

-0.11789 

-0.24663 

11 

-17.193 

-37.930 

0.08391 

0.00267 

12 

-11.397 

-  3 . 404 

0.04620 

0.01085 

13 

23.945 

42.668 

0.14752 

0.25374 

14 

55.349 

30.302 

0.27781 

0.23123 

15 

-10.846 

14.604 

-0.10033 

0.05556 

16 

59.557 

63.661 

0.37119 

0.14592 

17 

11.190 

2.355 

0.07027 

-0.03010 

18 

12.921 

-18.585 

0.09821 

0.04559 

19 

-17.522 

-10.792 

-0.29572 

0.02309 

20 

-39.395 

-19.506 

-0.13238 

-0.02975 

21 
22* 

-80.534 

-35.373 

-0.27386 

-0.03539 

23 

-12.884 

22.290 

-0.05130 

0.16258 

24 

-8.476 

51.982 

-0.03284 

0.27436 

25 

-25.229 

-12.730 

-0.13405 

-0.16915 

26 

22.956 

-0.910 

0.10552 

0.02366 

27 

-2.885 

-14.081 

-0.04947 

-0.13986 

28 

11.355 

7.095 

0.05162 

0.04954 

29 

30* 

-14.474 

2.247 

-0.07823 

0.00537 

31 

-44.763 

-21.355 

-0.19854 

-0.14849 

32 

-15.831 

25.363 

-0.08487 

0.08936 

33 

-11.938 

1.136 

-0.07265 

-0.00901 

34 
35* 

51.827 

15.167 

0.31837 

0.14705 

36 

-10.105 

-11.982 

-0.10263 

-0.06974 

37 

16.486 

8.212 

0.04258 

-0.01989 

38 

-6.060 

-41.126 

-0.00886 

-0.16374 

39 

9.809 

-32.049 

0.03889 

-0.00258 

40 

17.441 

-5.710 

0.13757 

-0.01437 

41 

50.106 

23.982 

0.19884 

0.01594 

42 

-15.341 

-22.285 

-0.02459 

-0.07543 

43 

44.717 

46.526 

0.26187 

0.26913 

44 

9.138 

28.242 

0.05684 

0.07427 

45 

-61.570 

-18.162 

-0.47302 

-0.20437 

46 

2.029 

-10.456 

0.02660 

-0.08615 

47 

-31.393 

-16.830 

-0.18230 

-0.17132 

48 

30.003  ' 

19.321 

0.24156 

0.15529 

*  Plots  eliminated  from  study  because  no  red  spruce  was  present. 


first  and  second  neighbors,  So  =  80,  S;l  =  138, 
and  S2  =  670. 

The  results  for  the  four  initial  residual  ring 
width  measures  are  summarized  in  table  11.  The 
spatial  autocorrelation  generally  appears  to  be 
higher,  based  on  the  second  nearest  neighbor 
weights.  Given  that  neighboring  plots  may 
sometimes  have  different  aspects  of  soil 
conditions,  it  can  be  assumed  the  set  of  weights 
based  on  the  first  and  second  order  nearest 
neighbors  gives  a  more  reliable  indication  of 
spatial  structure. 

Spatial  autocorrelation  coefficients  were  also 
calculated  for  the  residuals  from  the  regression 


analyses  shown  in  tables  9  and  10.  These  values 
are  presented  in  table  12.  The  standard  deviates 
were  computed  using  the  same  formulae  as  were 
used  to  obtain  table  11.  More  exact  expressions 
are  given  by  Cliff  and  Ord  (1981),  but  these  are 
very  tedious  to  use,  and  the  differences  in 
magnitude  are  generally  slight. 

Interpretation       of       Spatial       Analyses The 

results  in  table  11  show  clearly  that  pronounced 
spatial  dependence  exists,  whichever  measure  of 
residual  ring  width  is  used.  The  first  question 
that  arises  is  whether  such  dependence  can  be 
ascribed  to  purely  geographical  effects. 
However,   the   geographical   variables   do   not 
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Table  10 Residuals   from  weighted   least   squares   regression  for  RES20 , 

RES10,   PCT20,    and  PCT10 ,    respectively 


Observation 

RRES20 

RRES10 

PCRES20 

PCRES10 

1 

46.282 

-21.521 

0.24740 

-0.06286 

2 

-27.445 

-27.201 

-0.17263 

-0.16189 

3 

-35.316 

-22.902 

-0.16155 

-0.10963 

4 

5* 

6 

-14.947 

-45.974 

-0.28882 

-0.23249 

12.526 

12.103 

0.07637 

0.13770 

7 

-8.492 

-27.345 

-0.08494 

-0.12596 

8 

56.781 

77.765 

0.20111 

0.14621 

9 

-11.847 

-16.480 

-0.14614 

-0.14344 

10 

-19.367 

-33.018 

-0.16035 

-0.27380 

11 

-7.760 

-41.048 

0.11976 

-0.01109 

12 

-23.543 

-5.178 

0.06126 

0.02172 

13 

18.035 

38.016 

0.08884 

0.21764 

14 

57.645 

26.732 

0.25975 

0.22866 

15 

-18.072 

10.880 

-0.14416 

0.03250 

16 

37.223 

58.788 

0.19865 

0.09214 

17 

5.943 

-4.243 

0.05390 

-0.05542 

18 

23.991 

-23.934 

0.13074 

0.02145 

19 

-20.837 

-18.059 

-0.32424 

-0.00782 

20 

-42.430 

-26.236 

-0.15413 

-0.04206 

21 
22* 

-76.887 

-37.081 

-0.27356 

-0.03565 

23 

-17.231 

19.823 

-0.08004 

0.16041 

24 

-4.281 

51.511 

-0.03817 

0.25346 

25 

-37.022 

-14.344 

-0.13453 

-0.16752 

26 

28.356 

-0.332 

0.08779 

0.01210 

27 

-2.566 

-15.003 

-0.04214 

-0.15656 

28 

13.939 

6.884 

0.10000 

0.06398 

29 
30* 

-35.314 

-5.038 

-0.19339 

-0.01416 

31 

-48.186 

-28.191 

-0.20507 

-0.17886 

32 

-25.403 

26.303 

-0.09078 

0.12401 

33 

-23.449 

-6.225 

-0.13246 

-0.01701 

34 
35* 

56.338 

9.554 

0.34085 

0.12116 

36 

-18.279 

-14.585 

-0.12892 

-0.05920 

37 

17.310 

7.456 

0.02603 

-0.04263 

38 

-9.898 

-48.106 

-0.01409 

-0.16916 

39 

22.400 

-31.660 

0.09583 

0.02397 

40 

19.639 

-2.465 

0.12665 

-0.02988 

41 

44.136 

20.611 

0.15549 

-0.00942 

42 

-21.431 

-26.640 

-0.06461 

-0.10528 

43 

44.173 

43.168 

0.25985 

0.25656 

44 

5.930 

26.627 

0.04146 

0.05310 

45 

-59.933 

-15.349 

-0.47767 

-0.21879 

46 

-4.259 

-10.251 

0.00789 

-0.07634 

47 

-31.989 

-17.088 

-0.23592 

-0.16713 

48 

28.479 

19.917 

0.28237 

0.17932 

Plots  eliminated  from  study  because  no  red  spruce  was  presented. 


provide  any  degree  of  explanatory  power,  and  the 
spatial  pattern  of  the  residuals  is  essentially 
the  same.  The  next  question  is  whether  the 
biotic  factors  account  for  some  or  all  of  the 
spatial  structure. 

When  tables  11  and  12  are  compared,  it  may  be 
seen  that  the  level  of  spatial  autocorrelation 
has  diminished  in  all  cases.  For  the  variable 
PCT20,  the  autocorrelation  has  become  negative, 
but  this  may  be  due  to  the  uncorrected  effects  of 
autocorrelations  among  the  explanatory  variables. 
One  may  generally  conclude  that  the  regression 
analyses  account  for  much,  but  not  all,  of  the 
spatial  pattern  found  in  the  residual  ring  width 
values.    However,  one  should   recall  that  the 


major  variables  in  these  regressions  related  to 
stand  health.  While  this  analysis  indicates  that 
the  ring  width  changes  are  related  to  current 
stand  health,  the  reasons  for  that  current  health 
status  remain  undetermined. 

DISCUSSION  AND  SUMMARY 

The  results  of  steps  1  through  4  suggest  an 
overall  decrease  in  ring  width  increment  among 
red  spruce  in  the  Great  Smoky  Mountains  from  1965 
to  1975  or  1984,  relative  to  time  series 
forecasts  for  those  time  periods.  The  magnitude 
for  these  decreases  appeared  to  be  related  to 
elevation,  factors  of  stand  quality,  and   the 
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Table  11 . - -Results   of   tests   for  spatial   autocorrelation   among   original 
ring  width  residuals 

Variable"*" 


NN1      NN2  NN1     NN2 


RES20  0.256  0.425  1.48  3.21 

RES10  0.388  0.466  2.15  3.47 

PCT20  0.194  0.324  1.17  2.53 

PCT10  0.314  0.403  1.78  3.05 


+  RES20  and  RES10  refer  to  the  residuals  from  the  time  series  models  for 
ring  widths  averaged  over  the  forecast  periods,  20  and  10  years, 
respectively.  PCT20  and  PCT10  denote  RES20  and  RES10  divided  by  the 
overall  mean  ring  width  for  the  whole  series. 

NN1  and  NN2  refer  to  the  sets  of  weights  for  the  spatial  autocorrelation 
coefficient  based  upon  first  and  upon  first-  and  second- order  nearest 
neighbors,  respectively. 


Table  12 .- -Results   of   tests   for  spatial   autocorrelation   among  regression 
residuals**-' 


Variable          Coefficient               Standard  Deviate 
NN1     NN2 NN1 NN2 


Unweighted  Regression 


RES20 

0.018 

0.054 

PCT20 

-0.187 

-0.205 

RES  10 

0.277 

0.287 

PCT10 

0.093 

0.203 

Weighted 

Regression 

RES20 

-0.082 

-0.133 

PCT20 

-0.292 

-0.393 

RES10 

0.266 

0.268 

PCT10 

0.087 

0.202 

0.22  0.55 

■0.85  -1.28 

1.57  2.20 

0.60  1.59 

■0.30  -0.76 

■1.39  -2.59 

1.51  2.06 

0.57  1.57 


'See  table  11  for  definitions  of  variables,  coefficients,  and  standard 
deviates.   Regression  residuals  are  listed  in  tables  8  and  9. 
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relative  decreases  of  first  and  second  order 
nearest  neighboring  plots.  These  findings  should 
be  compared  to  those  generated  by  other 
approaches  to  gain  an  understanding  of  the  impact 
of  certain  key  decisions  in  the  stages  of  the 
analysis . 

Both  Landau  and  others  (1985)  and  Cook  (1987) 
recommended  that  annual  basal  area  growth 
increment  is  preferable  to  ring  width  as  a 
measure  of  annual  productivity.  Use  of  this 
measure  could  have  reduced  the  nonstationarity  in 
several  of  the  time  series .  One  approach  to  time 
series  that  display  marked  trends  is  to  transform 
the  data  (Cook  1987).  Instead,  we  followed  the 
usual  ARIMA  paradigm  and  differenced  the  series 
where  there  were  marked  trends ;  20  of  the  44 
series  actually  required  differencing.  It  would 
be  useful  to  compare  forecasts  obtained  from  the 
untrans formed,  but  possibly  differenced,  series 
to  those  obtained  from  transformed  series . 
Additionally,  the  average  of  the  ring  width  of 
the  trees'  cores  that  met  our  inclusion  criteria 
was  used.  Landau  and  others  (1985)  recommended 
the  use  of  trimmed  means.  Specific  circumstances 
of  other  data  sets  may  determine  the  advisability 
of  one  procedure  over  another  as  the  best  way  to 
minimize  heterogeneity. 

The  methods  used  in  each  of  these  steps  are 
capable  of  further  refinement,  and  several 
suggestions  are  included  in  the  Recommendations 
section.  Nevertheless,  the  basic  paradigm 
represents  a  substantive  approach  to  the 
evaluation  of  recent  trends  in  the  width  of  tree 
rings . 

Results  obtained  from  the  steps  taken  in  this 
study  should  be  compared  with  results  generated 
from  other  approaches.  Future  studies  might 
include  the  use  of  basal  area  increment  rather 
than  ring  width  as  a  measure  of  annual 
productivity  and  may  incorporate  transfer 
function  models  of  climatic  factors  as  well  as 
intervention  analysis  to  filter  out  the  effect  of 
important  forest  perturbations. 

RECOMMENDATIONS 

1.  Following  Landau  and  others  (1985),  it  is 
recommended  that  future  studies  should  use 
incremental  basal  area  rather  than  ring 
width. 

2.  The  possibility  of  using  trimmed  means  rather 
than  arithmetic  means  should  be  considered. 
However,  the  issue  of  how  well  measurements 
on  five  healthy  trees  reflect  overall  stand 
health  requires  further  examination. 

3.  It  is  important  to  look  for  changes  in  ring 
width  or  other  indicators  relative  to  what 
might  be  expected.  The  forecasting  approach 
used  in  this  report  is  one  way  of  excluding 
such  trends,  but  others  should  be  examined. 
The  study  of  proportional  changes  seems 
preferable  to  that  of  absolute  changes. 

4.  In  future  time  series  analyses,  automated 
procedures  might  be  used  (AUT0B0X,  developed 
by  David  Reilly  of  Automatic  Forecasting 
Systems ,  Inc . ) . 

5 .  Where  known  problems  of  fires ,  aphids ,  or 
infestations  occur  on  particular  plots, 
intervention  analysis  should  be  used  to 
filter  out  such  effects. 


6.  The  inclusion  of  biotic  variables  in  the 
regression  models  serves  to  link  the  change 
indicators  to  stand  health.  However,  it  does 
not  resolve  how  the  stands  came  to  be  in  that 
condition.  The  lack  of  any  worthwhile 
correlations  between  the  indicators  and  the 
locational  variables  suggests  that  other 
factors  may  be  at  work  in  determining  stand 
health;  furthermore,  the  high  levels  of 
spatial  autocorrelation  in  the  ring  width 
change  data  indicated  that  such  factors  are 
spatially  concentrated. 

7.  The  extremely  variable  topography  and 
locations  of  the  sites  suggest  that  purely 
spatial  models  (Cliff  and  Ord  1981)  are 
unlikely  to  be  of  direct  value  in  this 
study.  However,  the  potential  exists  for 
worthwhile  applications  with  more  homogeneous 
clusters  of  sites. 

8.  Future   analyses   could   include   transfer 
function   models   involving   climatic 
variables,   once  detailed  models  of  these 
variables  have  been  developed. 
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APPENDIX 


FORTRAN  program  for  computing  the  spatial 
autocorrelation  coefficient 


(A)   Program 

DIMENSION  Z(100),  X(100) ,  NPLOT(IOO) ,  NNl(lOO) 
NN2(100),  *WT1(100),  WT2(100,  NCROSS(250) 

READ,  N 

READ,  SOI,  Sll,  S21,  S02 ,  S12,  S22 

SUMX-O . 0 

DO  10  I-l.N 

READ,  NP,  Nl,  N2,  NREF 

NPLOT  (I)-NP 

NN1(I)-N1 

NN2(I)-N2 
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NCROSS(NP)-=NREF 
10   CONTINUE 

DO  15  I-l.N 

READ,  XA 

SUMX=SUMX+XA 

X(I)=XA 
15    CONTINUE 

SUMZZ=SUMZ4=0 . 0 

SAC1=SAC2=0.0 

XBAR=SUMX/N 

DO  20  1=1,  N 

Z(I)=X(I)-XBAR 

SUMZZ=SUMZZ+Z ( I ) **2 

SUMZ4=SUMZ4+Z ( I ) **4 
20     CONTINUE 

B2=N*SUMZ4/ ( SUMZZ**2 ) 

PRINT,  'KURTOSIS=' ,B2 

DO  30  J=1,N 

KA=NN1(J) 

KB=NN2(J) 

KC=NPLOT(J) 

IF   (KA.EQ.O)  GO  TO  30 

JA=NCROSS(KA) 

JC=NCROSS(KC) 

SAC1=SAC1+Z(JC)*Z(JA) 

IF  (KB.EQ.O)  GO  TO  30 

JB=NCROSS(KB) 

SAC2=SAC2+Z(JC)*Z(JB) 
30     CONTINUE 

SAC2=SAC1+SAC2 

SAC1=SAC1/SUMZZ 

SAC2=SAC2/SUMZZ 

SD1=N*((N*N-3*N+3)*S11-N*S21+3*S01*S01) 

SD1=SDI-B2*((N*N-N)*S11-2*N*S21+6*S01*S01) 

SD1=SD1/((N-1)*(N-2)*(N-3)*S01*S01) 

SD1=SD1-1.0/(N-1)**2 

SD1=SQRT(SD1) 

SD2=N*((N*N-3*N+3)*S12-N*S22+3*S02*S02) 

SD2=SD2-B2*((N*N-N)*S12-2*N*S22  +6*S02*S02) 

SD2=SD2/((N-1)*(N-2)*(N-3)*S02*S02) 

SD2=SD2-1.0/(N-1)**2 

SD2=SQRT(SD2) 


PRINT, 'SD1=' ,SD1, 'SD2=' , SD2 

PRINT,'  SPATIAL  A/C  FOR  FIRST  NN  IS 

PRINT,'  SPATIAL  A/C  FOR  SECOND  NN  IS 

SAC1=(SAC1+1 . 0/(N-l) )/SDl 

SAC2=(SAC2+1  0/(N-l))/SD2 

PRINT, 'STD.  DEVIATE  FOR  FIRST  NN  IS  ' 

PRINT, 'STD.  DEVIATE  FOR  SECOND  NN  IS 

STOP 

END 


,SAC1 
' , SAC2 


,SAC1 
' , SAC2 


(B)   Plot  Number- -first  nearest  neighbor- -second 
nearest  neighbor--  order  of  plot  in  listing  of 
values  (required  inputs  to  vectors  NPLOT  NN1 
NN2,  AND  NCROSS) 


32  54  55  29 

54  32  55  35 

55  54  32  36 
7  159  207  6 


7  43 

7  5 

159  44  125  4  0 

159  41  4  125  0 


39 
4 


159  6 
6  159 
207  7 
153  6 
12  9  25  11 

9  25  12  8 

25  9  12  23 

26  27  25  24 

27  26  10  25  46  17  8  34 

10  11  113  9 

113  10  11  38  8  17  46 

11  10  113  10  19  18  0 
21  11  10  20 

28  23  11  26 

23  24  28  21 

24  23  14  22 
14  24  23  13 
16  109  14  15 
109  16  14  37 


7 
18 


(continued) 
15  13  37  14 
158  0  0  42 
13  37  36  12 
37  36  13  33 
36  37  13  32 
132  0  0  40 


34  29  33  31 
29  33  34  27 
33  29  34  30 
31  46  8  28 

17  8  46  16 


18  19  0  17 

2  20  1  2 
20  1  2  19 
1  20  2  1 

3  1  20  3 
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A  Fractal  Approach  to  Analysis  of  Tree  Ring  Increments 

R.  A.  J.  Taylor 


SUMMARY 

Information  from  plots  established  in  the  Great 
Smoky  Mountains  by  the  National  Park  Service  and 
the  Tennessee  Valley  Authority  was  used  to 
determine  annual  tree  ring  widths  from  core 
samples  of  red  spruce  (Picea  Rubens  Sarg.).  The 
red  spruce  core  samples  showed  a  significant 
dependence  of  variance  on  the  mean  size  of  tree 
rings  at  67  of  68  plots.  At  9  of  48  plots,  the 
variance  has  increased  more  rapidly  since  1943; 
of  the  others,  7  have  shown  a  decrease  since  1940 
and  32  showed  no  change.  The  dependence  of 
variance  on  mean  of  a  measurement  was  interpreted 
in  terms  of  "fractals,"  a  term  coined  to  denote 
fractional  dimension.  The  change  in  fractional 
dimension  over  time  indicated  an  evolution  of 
factors  that  influenced  the  dependence  of 
variance  on  mean;  these  factors  may  have  been 
successional ,  climatic,  or  anthropogenic,  all  of 
which  seemed  to  vary  on  about  the  same  time 
scale.  It  was  concluded  that  variance -mean 
analysis  may  be  an  inexpensive  and  promising  area 
of  inquiry  in  dendrochronology. 


INTRODUCTION 

Mortality  of  large  forest  areas  in  several 
parts  of  the  world  has  caused  fear  that  the 
concentration  of  anthropogenic  compounds  in  the 
atmosphere  may  be  increasing.  Acids  and  other 
oxidizing  agents  of  human  origin,  notably  ozone, 
have  been  detected  in  the  atmosphere  and  are 
probably  capable  of  interfering  with  tree  growth 
and  survival.  However,  there  are  no  data 
concerning  the  level  of  atmospheric  pollution 
(Cook  1987,  Kiester  and  others  1985).  Therefore 
other  possible  influences  and  causes  must  be 
considered,  such  as  the  effect  of  climate  on  tree 
ring  growth.  Because  climate  and  anthropogenic 
effects  are  easily  confounded,  this  report  will 
focus  on  the  analysis  of  change  and  not  on 
distinguishing  between  pollution  and  climate. 

Annual  tree  rings  in  temperate  regions  provide 
a  convenient  record  of  a  tree's  growth  history. 
Comparison  of  tree  rings  over  a  geographical  area 
has  frequently  been  used  to  determine  climatic 
changes;  it  is  assumed  that  patterns  common  to 
all  trees  of  the  same  species,  similar  age,  and 
in  the  same  soils  should  respond  alike  to  weather 
conditions  that  are  basically  unvarying  (Creber 
1977,  Guiot  and  others  1982). 

BACKGROUND 

It  is  commonly  assumed  that  the  width,  W(t),  of 
a  tree  ring  laid  down  in  year  t  is  the  linear  sum 
of  four  systematic  components  and  a  random 
component: 


W(t)  =  A(t)  +  B(t)  +  C(t)  +  D(t)  +  e(t) 
where 

A(t)  =  age  factor (s)  unique  to  each  tree, 

B(t)  =  disturbances  unique  to  each  tree, 

C(t)  =  climatic  effects  common  to  all  trees  at 

a  site, 
D(t)  =  disturbances  common  to  all  trees  at  a 

site,  and 
e(t)  =  random  component. 

It  may  be  helpful  to  review  some  of  the 
assumptions  frequently  made  about  W(t) .  The 
successive  increments  are  assumed  to  be 
identically  distributed.  For  the  benefit  of 
certain  analyses ,  the  increments  are  further 
assumed  to  be  statistically  independent,  with  the 
marginal  distribution,  e(t),  Gaussian  with  zero 
mean,  and  constant  variance.  Such  a  time  series 
is  called  a  stationary  Gaussian  random  walk  or 
Brownian  motion.  While  no  one  working  in 
dendrochronology  seriously  expects  the  Brownian 
assumption  to  be  valid,  the  nature  of  statistics 
often  demands  that  it  be  assumed. 

The  concepts  of  randomness  are  context 
dependent,  and  the  word  may  be  used  in  a 
confusing  manner.  There  are  two  broad, 
alternative  definitions  of  randomness: 
"predictable  behavior,  efficiently  described  by  a 
statistical  probability  distribution"  and 
"haphazard  behavior,  governed  by  no  known  rules" 
(Mandelbrot  1967).  In  addition,  events  in  which 
the  mean  and  variance  are  equal  (Poisson 
distributed)  are  often  said  to  be  random. 
Predictable  behavior  is  synonymous  with 
stochastic  and  differs  from  deterministic  because 
its  expected  or  average  outcome  is  predictable, 
while  the  specific  outcome  of  a  trial  lies  within 
certain  bounds  defined  by  the  variance. 

As  stated  above,  statistical  independence  of 
successive  increments  is  a  well-known 
simplification.  The  assumption  of  stationarity , 
however,  has  a  special  implication  that  is  rarely 
questioned:  the  sample  moments  vary  little  from 
sample  to  sample ,  provided  the  samples  are  large 
enough.  In  our  analysis,  this  assumption  is 
relaxed,  and  the  assumption  is  made  that  the 
variance  is  infinite  or  at  least  so  large  that  it 
may  be  treated  as  infinite.  Thus  the  assumption 
of  Gaussian  marginal  distribution  is  abandoned  in 
favor  of  the  Cauchy,  permitting  the  Central  Limit 
Theorem  to  be  invoked  without  assuming  constant 
variance  (Mandelbrot  1969) .  Therefore  the  need 
to  amend  the  other  assumptions  of  independence 
and  stationarity  is  alleviated  (Berger  and 
Mandelbrot  1963) . 

The  assumption  of  infinite  variance  is 
equivalent  to  assuming  randomness  that  is 
predictable  but  haphazard;  the  long-term  trend  is 
evident,  but  the  short-term  signal  is 
unpredictable.    Theory  expounded  by  Mandelbrot 
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(1960,  1963,  1969)  related  haphazard  time  series 
with  power  laws  and  distributions  with  infinite 
variance . 

THEORETICAL  METHODS 

Consider  the  series  x(t)  ,  t  =  1,2 n,  of 

observations  taken  at  equal  intervals  of  time  or 
space.  Any  particular  series  [x(t)]  is  assumed 
to  be  the  realization  of  a  process  W(t)  that  will 
be  defined  later;  t  is  used  as  an  indicator 
variable  at  equal  intervals  of  time  or  space. 
Defining  the  following  variables: 

V(k)  -  E[x(t)-x(t+k;>]2  -  2[V(x)-C(x;k)]  , 

V(x)  =  E[x(t)-M(x)]2, 

M(x)  =  E[x(t)], 

C(x;k)  =  E{[x(t)-M(x)][x(t+k)-M(x)]}. 

V(x)  is  the  variance  of  the  series  with  M(x) 
the  mean.  Theory  based  on  the  usual 
interpretation  of  the  Central  Limit  Theorem 
permits  one  to  assume  that  these  parameters 
estimate  population  values.  However,  in  this 
analysis  the  assumption  is  not  necessarily  valid: 
the  parameters  simply  represent  the  sample  values 
and  are  therefore  not  asymptotic  to  the 
population  values.  C(x;k)  is  the 
spatial/temporal  covariance  across  the  data  and 
is  related  to  the  standard  treatments  of 
time/space  series  data,  including  tree  ring 
analyses  (Guiot  and  others  1982) .  Information  is 
given  in  C(x;k)  on  the  regular  variations  in  the 
data  at  periods  equal  to  k.  The  variogram  is 
V(k)  and  provides  information  on  the  nonregular 
variation  at  lag  k. 

Consider  V(k)  as  k  varies:  dividing  V(k)  = 
2{V(x)-C(x;k)}  through  by  V(x)  gives  V(k)  = 
2V(x){l-R(k) } ,  where  R(k)  is  the  serial 
correlation  coefficient  that  takes  values 
-l<R(k)<l;  therefore  V(k)=2rV(x) {1-r } ,  0<r<2 ,  and 
0<V(k)<V(x).   Thus: 

V(k)  =  0  when  serial  correlation  at  lag  k  is 

1, 
V(k)  =  2V(x)  when  no  serial  correlation  at  lag 

k,  and 
V(k)  =  4V(x)  when  serial  correlation  at  lag  k 

is  -1. 
At  any  specific  value  of  k,  say  k*,  V(k)  will 
give  information  on  all  variations  not  having  a 
cycle  at  k. 

To  obtain  information  on  all  variations,  V(x) 
is  commonly  used,  but  we  cannot  say  what 
variation  we  have  at  t  =  i,  only  that  it  is  over 

the  interval  t  =  i ,  i+1 ,  i+2 and  so  forth. 

Ideally,  the  variance  would  be  partitioned  in  the 
manner  of  V(k)  but  without  limitation.  If  k  goes 
to  zero,  all  scales  of  variation  are  included  and 
simultaneously  V(k)  disappears.  To  estimate 
V(0)  ,  compute  V(k)  for  k>0  and  extrapolate  back 
to  zero.  To  obtain  the  rate  of  approach  of  V(k) 
to  the  origin,  plot  logV(k)  against  log(k)  and 
determine  the  gradient,  dlogV(k)/dlog(k)  =  ft  as 
k— >0 . 

Another  approach  to  estimate  the  gradient  of 
V(0)  is  to  replicate  the  generating  process  W(t): 

Wi(t),  W2<t) ;  V(0)  can  then  be  estimated  at 

any  or  all  t.  The  estimates  V^(0)  are  all  at  the 
origin,  so  the  gradient  must  be  extracted  from 
them.  Plotting  the  means  and  variances  of 
several  series  against  time  shows  that  the  actual 


values  vary  according  to  no  particular  pattern; 
however,  it  must  be  noted  that  the  magnitude  of 
Vt(0)  seems  to  increase  with  time  (figs.  1-8). 
Since  Vt(0)  is  partly  dependent  on  Mt(0) , 
standardizing  Vt(0)  with  Mt(0)  may  show  a  trend. 
Figures  1  through  8  also  show  the  coefficient  of 
variation  (yvt(0)/Mt(0) )  increasing  with  time. 

Figure  9  plots  logV(t)  against  logM(t)  and 
shows  how  the  variance  increases  with  respect  to 
mean.  The  gradient  dlogV/dlogM  -  b  is  an 
estimate  of  ft.  This  is  easily  demonstrated  by 
the  following  argument.  Define  a  reference  mean 
Mq  and  a  comparison  mean  Ms  =  sMq  (s>1} .  Now  Vq 
=  aM0b  and  Vs  =  aMsb  =  a(sM0)b  =  sb(am0b)  =  sbV0. 
Evidently  the  variance  is  scale  independent,  and 
its  gradient  b  is  an  intrinsic  component,  as 
expected  from  dlogV(k)/dlog(k) . 

Taking  logs:  logVs  =  logVQ  +  blog(s) ,  remember 
that  Vq  is  a  variance  corresponding  to  an 
arbitrarily  chosen  mean  Mq.  Vq  is  thus  also 
arbitrary;  Mq  can  be  chosen  such  that  Vq  is  1. 
Thus  logVs  =  blog(s) ,  which  is  very  nearly 
logV(k)  =  £log(k) .  Although  k  and  s  are  not  quite 
synonymous,  if  multiples  are  chosen  as  values  of 
s,  then  logV  =  blog(i)  +  C  describes  the  locus  of 
variance  at  spatial/temporal  intervals  i  = 
1,2,...,  and  is,  except  for'  C,  identical 
logV(k)  =  ftlogK.   Therefore  &  =  b. 

The  same  result  can  be  obtained  from 
empirical  argument.  Let  there  be  a  series 
samples  taken  along  a  transect  AB.  Divide 
into  NK  intervals  such  that  there  are  N  groups  of 
K  intervals.  Neither  the  K  groups  nor  the  N 
groups  need  to  be  contiguous,  but  for  simplicity 
it  is  assumed  that  they  are.  The  number  of 
entities  in  each  of  the  NK  intervals  is  counted. 
In  the  present  case,  the  size  of  the  tree  ring 
increment  is  measured:   xj*  where  x  is  the  size 

and  i  =  1 , 2 , . . . , N  and  j  =  1 , 2 K.   Now  compute 

the  mean  Mi  and  variance  V^  for  each  of  the  i  = 
1,2, ... ,N  groups  from 


to 

an 
of 
AB 


Mj.  =  Zxy/K 


Vi  =  S(xij-Mi)V(K-l) 


Let  the  central  interval  of  each  group  become 
the  center  of  mass  for  that  group.  The  center  of 
mass  has  at  least  two  parameters  describing  the 
distribution  of  observations  within  the  group: 
Mt  and  V^ 

If  logV  is  plotted  against  logM,  empirically 
they  are  related  by  the  power  law,  V  =  aM  ,  where 
a  and  b  are  empirically  estimable  parameters 
(Taylor  1961,  Perry  1981).  Consider  the 
arbitrary  series  x,  (j  =  1,2,...,K)  that  has 
mean  and  variance  M  and  V  ;  now  compute  the 
serial  covariance  C  (k)  having  lag  k,  from 

C*(k)  =  E^xj*  -  M*][xj+k*  -  M*]. 
From  V   and  C  (k)  the  variance  of  increments  is 
computed: 


xj+k*)2  -  2[V* 


C*(k) 


V*(k)  -  E(xj* 

Half  V  (k)  is  referred  to  as  the  variogram,  and 
its  computation  is  the  first  step  in  the 
interpolation  process  known  as  kriging  (Journel 
and  Huijbregts  1978).  The  covariance  term  is 
C  (k)  and  is  related  to  the  systematic  part  of 
the  ring  increments,  particularly  growth  at  small 
k.  At  larger  k,  long-period  variations  such  as 
climatic  cycles  predominate.  Filtering 
techniques  (Guiot  and  others  1982)  concentrate  on 
the  structure  of  C  (k)  over  varying  periods  of 
time.  The  total  variance  of  the  series  V 
contains  both  the  systematic  and  nonsystematic 
variance. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  1  (1850  —  1984) 
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Figure  l.--rime  series    of    mean,     variance ,     and    coefficient    of    variation    of 
tree-ring  increments  from  1850   to  1984  at  Tennessee  Valley 
Authority  site  1. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  9  (1850  —  1984) 
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Figure  2. --Time   series   of  mean,    variance,    and  coefficient   of  variation   of 

tree-ring     increments     from     1850      to     198U     at     Tennessee     Valley 
Authority  site  9. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  23  (1850  —  1984) 
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Figure  3. --Time  series  of  mean,    variance ,    and  coefficient  of  variation  of 
tree-ring   increments   from  1850    to   1984   at   Tennessee   Valley 
Authority  site  23. 
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MEAN  ANNUAL  TREE-RING  WIDTH  -  SITE  36  (1850  -  1984) 
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Figure  U.--Time  series  of  mean,    variance ,    and  coefficient  of  variation  of 
tree-ring  increments  from  1850   to  1984  at  Tennessee  Valley 
Authority  site  36. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  307  (1901  —  1984) 
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Figure  5. --Time  series  of  mean,    variance,    and  coefficient  of  variation  of 

tree-ring    increments    from    1850    to    1984    at    National    Park    Service 
site  307. 


46 


MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  310  (1901  —  1984) 
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Figure  6. --Time  series  of  mean,    variance ,    and  coefficient  of  variation  of 
tree-ring  increments  from  1850   to  1984  at  National  Park  Service 
site  310. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  316  (1901  —  1984) 
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Figure  7. --Time    series    of    mean,     variance ,     and    coefficient    of    variation    of 
tree-ring  increments  from  1850   to  19&U  at  National  Park  Service 
site  316. 
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MEAN  ANNUAL  TREE— RING  WIDTH  —  SITE  321  (1901  —  1984) 
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Figure  8. --Time  series  of  mean,    variance,    and  coefficient  of  variation  of 
tree-ring  increments  from  1850   to  198U  at  National  Park  Service 
site  321. 
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This  study  focuses  primarily  on  the  systematic 
variance  common  to  all  trees  at  a  site. 
Therefore  the  variance  across  all  cores  at  a  site 
in  each  year  was  calculated  to  form  the  series 
Vt.  The  common  systematic  components  at  a  site 
formed  a  baseline  from  which  all  other  components 
of  variance  were  referenced.  To  find  the  common 
systematic  component,  variance  was  plotted 
against  mean  to  standardize  variance.  The  mean 
was  assumed  to  be  linearly  related  to  the 
baseline. 

Also  to  be  noted  is  the  change  in  the  common 
component,  represented  by  the  change  in  Vt 
relative  to  Mt.  Changes  between  the  relationship 
of  Vt  and  Mt  indicate  changes  in  the  normal 
behavior  of  the  process  W(t).  Changes  result 
from  evolution  in  the  process  variance.  If  a 
long-term  change  is  anticipated,  then  evolution 
in  the  system  process  can  be  demonstrated  by 
determining  differences  in  the  rates  of  change  of 
variance  before  and  after  the  reference  point. 


PROCEDURE 

Of  the  plots  established  by  the  National  Park 
Service  (NPS)  and  the  Tennessee  Valley  Authority 
(TVA)  ,  only  those  with  five  or  more  red  spruce 
trees  in  the  sample  were  selected  for  analysis. 
With  2  cores  (series)  per  tree,  means  and 
variances  for  each  year  from  the  beginning  of  the 
series  were  computed  from  10  ring  width 
estimates.  Series  ranged  in  length  from  40  to 
135  years.  These  data  yielded  20  bivariate 
series  with  length  of  approximately  84  years  (NPS 
data) ,  and  48  series  with  lengths  varying  from  40 
to  135  years  (TVA  data) .  The  logarithm  of 
variance  was  then  regressed  on  log  mean  with  each 
point  representing  1  year. 

RESULTS 

Figures  1  through  8  show  the  mean  (M) ,  variance 
(V) ,  and  coefficient  of  variation  (CV)  against 
time  at  TVA  sites  1,  9,  23,  and  36  and  NPS  sites 
307,  310,  316,  and  321.  Both  the  mean  and 
variance  were  highly  variable,  and  the  plots  of 
CV  against  time  gave  the  strong  impression  that 
the  variance  was  increasing  with  time. 
Furthermore,  the  detailed  behavior  differed 
greatly  from  site  to  site.  Figures  9  through  16 
show  variance -mean  regression  plots  for  the  same 
8  sites,  and  table  1  shows  the  results  of 
regressions  of  68  sites.  Sixty-seven  regressions 
were  significant  at  probabilities  of  less  than 
0.05.  Only  site  113  showed  no  dependence  of 
variance  on  mean. 

To  investigate  the  possibility  that  conditions 
in  the  post-World  War  II  era  were  different  from 
prewar  conditions,  the  data  of  sites  with  runs 
longer  than  80  years  were  split  at  1943  and  the 
variance-mean  regressions  repeated  to  test  the 
hypothesis  that  the  acceleration  of  variance  with 
mean  has  changed.  Table  2  shows  that  of  the  28 
TVA  data  sets  suitable  for  analysis,  8  showed 
increases  in  dependence,  18  showed  no  change,  and 
none  showed  a  decrease  in  dependence.  There  were 
two  sites  with  a  nonsignificant  regression. 
Interestingly,  the  NPS  data  showed  the  reverse 
pattern:  only  one  site  showed  a  significant 
increase  in  variance-mean  dependence  since  1943, 


nine  sites  showed  no  difference  and  seven  showed 
a  significant  decrease  since  1943,  three  sites 
had  a  nonsignificant  regression. 

Preliminary  analyses  with  multivariate  methods 
(principal  components  analysis ,  PCA)  failed  to 
identify  any  differences  in  regression  gradient 
attributable  to  differences  in  the  two  data  sets. 
Variables  included  in  the  PCA  were  annual 
rainfall,  mean  annual  maximum  and  minimum 
temperatures,  xeric/mesic  status,  altitude,  and 
stand  basal  area. 

DISCUSSION 

Tree  rings  reflect  the  net  effect  of  age, 
health,  soil,  and  biotic  conditions  on  tree 
growth . 

Age  and  health  are  unique  aspects  of  each 
individual  tree,  but  the  environmental 
conditions  experienced  by  all  trees  in  close 
proximity  are  considered  similar  or  common 
external  factors.  However,  competition  is  one 
aspect  of  a  tree's  external  environment  that  is 
considered  unique. 

Therefore  three  influences  of  growth  can  be 
established: 

1.  Age/growth — unique  internal  factors, 

2.  Competition — unique  external  factors,  and 

3.  Environment — common  external  factors. 

When  samples  are  composed  of  mature  and/or  old 
trees,  the  age - dependent  variance  will  be 
relatively  small  and  may  even  be  missing  in  the 
younger  trees .  When  very  young  and  very  old 
trees  are  included  in  the  sample,  age/growth  is 
likely  to  have  an  impact  only  in  the  early  part 
of  a  series  and  variances  are  higher  than  normal. 

Competition  between  trees  in  the  sample  may 
help  one  to  understand  the  dependence  of  variance 
on  mean.  For  a  given  value  of  the  mean,  a  high 
variance  indicates  a  wider  range  of  individual 
responses .  As  the  mean  increases ,  the  dependence 
of  variance  on  the  mean  suggests  that  only  some 
elements  of  a  sample  are  increasing,  while  others 
may  be  decreasing  or  increasing  very  little.  If 
increases  in  tree  ring  width  were  all  equal  at  a 
site,  the  variance  would  be  increasing 
proportionately  and  b  -  1.  The  mean  value  of  b 
was  1.77,  suggesting  that  at  the  majority  of 
sites,  when  conditions  were  conducive  to  growth 
increases,  only  some  trees  responded.  A  change 
in  b  over  time  was  interpreted  as  a  change  in  the 
relative  competitiveness  of  the  trees  at  a  site 
due  to  a  change  in  the  external  conditions. 

CONCLUSIONS 

Because  the  relationship  between  variance  and 
mean  is  still  uncertain,  no  definite  conclusions 
can  be  made  about  the  influences  on  tree  growth. 
Because  the  nature  of  fractional  dimension  is 
still  poorly  understood,  the  theoretical 
treatment  is  not  rigorous.  Only  since  the 
publication  of  Mandelbrot's  book  (1982)  has  there 
been  an  increased  awareness  and  interest  in 
fractals  among  scientists  other  than  topologists. 

The  ideas  developed  in  this  analysis  were 
preliminary,  but  the  reduction  of  the  masses  of 
tree  core  data  to  a  single  parameter,  b,  was  a 
new  contribution  to  statistical  methods  in 
dendrochronology.  However,  the  results  of 
splitting   the   data   into   two   sections   were 
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Figure  9 .- -Variance-mean  plot  (log  scale)  of 
Tennessee  Valley  Authority  site  1. 
All  regressions  are  significant  at 
p<0.01.  This       plot       shows       some 

nonlinear ity. 


Figure  10 .--Variance-mean   plot    (log  scale)   of 
Tennessee   Valley  Authority   site      9. 
All    regressions    are    significant    at 
p<0.01.         This    plot    shows    very    good 
regression . 


SITE  23 


SITE  36 


600 
4.75 
4.50 
425 
4.00 
3.75 
3.50 
3.25 
300 
2.75 
2.50 
US 
200 


o 

oo    o         o 

oo              > 

o  o  oo  y 

o 

ooo  >^b  o 

oo 

o 

s€       00 

o  00     o 

ro       o  o      < 

o  . 

o    o 

o     o>o 

o 

ooo 

o    o 

y^  o 

o 

0 

*  o 

o 

oo 

o 

0 

0  o 

o  oo  oo 
00  o 


1.7    1.6    1.9    20    2.1    2.2    2.3 
LOG (MEAN) 


2.4 


25 


4.7S 

4.50 

4.2S 

400 

bJ 
o 

z 
< 

375 
350 

_l 
26 


£  325 

> 

O  300 

o 

ml 

2.75  - 
2.50 
£25  - 

200  - 
1.75 


o                    / 

o 

ooo           of 

oo  0    db 

oo  o  o    ooo  y  00 

o 

o  o     o    y  o 

o  o  oo       y            0  0 

o       oo  jr     o    oo  o 

O       000  OyOO     O 

O          O  or                     o 

o 

O        O       /        O       O 

/          00  O    O 

o  o    oyb  o         oo 

Of  OO               o 

Or  oo    o      o 

18 


L9 


_l_ 


20 


2.1      2.2 
LOG (MEAN) 


2.3 


2.4 


8.5 


Figure  11 .- -Variance -mean  plot  (log  scale)  of 
Tennessee  Valley  Authority  site  23. 
All  regressions  are  significant  at 
p<0.01.  This  plot  shows  very  good 
regression. 


Figure  12 .- -Variance-mean  plot  (log  scale)  of 
Tennessee  Valley  Authority  site  36. 
All  regressions  are  significant  at 
p<0.01.  This  plot  shows  very  good 
r eg res sion.  51 
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Figure  13 .- -Variance-mean  plots  (log  scale)  of 
National  Park  Service  site  307.  All 
regressions  are  significant  at  p<0.01. 
This  plot  shows  very  good  regression. 


Figure  14 .- -Variance-mean  plots  (log  scale)  of 
National  Park  Service  site  310.  All 
regressions  are  significant  at  p<0.01. 
This  plot  shows  very  good  regression. 
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Figure  15 .- -Variance-mean  plots  (log-scale)  of 
National  Park  Service  site  316.  All 
regressions  are  significant  at  p<0.01. 
This  plot  shows   some  nonlinear ity . 


Figure  16 .- -Variance-mean  plots  (log  scale)  of 
National  Park  Service  site  321.  All 
regressions  are  significant  at  p<0.01. 
This  plot  shows  very  good  regression. 
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inconclusive.  A  reduction  in  variance  in  recent 
years  at  each  TVA  site  may  be  indicated  by  the 
TVA  data,  while  the  NPS  sites  have  become  more 
variable.  If  this  is  so,  the  difference  between 
the  two  data  sets  must  be  determined.  The 
results  of  the  PCA  were  also  inconclusive,  but 
the  set  of  site  characteristics  used  in  the 
analysis  was  very  small. 


The  fractal  approach  to  the  analysis  of  tree 
ring  widths  is  a  promising  area  for  further 
research.  However,  this  method  may  not  help  to 
identify  anthropogenic  influences  on  tree  growth; 
change  may  be  determined,  but  the  cause  of  that 
change  may  still  require  carefully  controlled, 
long-term,  large-scale  forest  experiments . 


Table  1 . - -Regression   analyses   of   Che   log    (variance)    against    log    (mean)    of 
tree  ring  increments* 


Site 
No.    N 


Log(a) 


Standard   Adjusted 
error  of  b     R2 


F- 
ratio 


Significance 


Tennessee  Valley  Authority  Plots 


1 

135 

1.38 

1.05 

2 

71 

-0.72 

1.84 

3 

88 

1.50 

0.96 

4 

92 

-1.98 

2.55 

5 

35 

-1.04 

1.90 

6 

96 

1.83 

0.77 

7 
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0.41 

1.55 

8 

57 

-1.02 

2.05 

9 
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-0.32 

1.85 

10 
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1.18 

1.06 

11 

45 
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0.94 

12 
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1.64 

13 

45 

-3.22 

3.08 

14 
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15 
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Table  1 .- -Regression   analyses   of   the   log    (variance)    against    log    (mean)    of 
tree  ring  increments- -Continued. 


Site 
No.    N 


Log (a) 


Standard   Adjusted 
error  of  b 


R/ 


F- 

ratio 


Significance 


National  Parks  Service  Plots 


301 

84 

-1.50 

2.60 

0.187 

0.70 

193.0 

302 

84 

-0.88 

2.05 

0.205 

0.54 

100.0 

303 

84 

-1.57 

2.58 

0.201 

0.66 

164.0 

304 

84 

2.04 

0.75 

0.237 

0.10 

10.1 

** 

305 

84 

0.04 

1.63 

0.111 

0.72 

214.0 

307 

84 

-0.20 

1.88 

0.121 

0.74 

242.0 

308 

84 

-1.75 

2.62 

0.352 

0.40 

55.4 

309 
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-2.46 

2.88 

0.332 

0.47 

75.1 
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84 

-1.25 

2.38 

0.150 

0.75 

253.0 
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84 

0.05 

1.66 
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0.55 

103.0 
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84 
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** 
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2.54 

0.207 

0.64 

152.0 

316 

84 

1.63 

0.73 

0.244 

0.10 

0.9 

** 

317 

84 

-0.13 

1.72 

0.169 

0.55 

104.0 

318 

84 

-1.23 

2.36 

0.250 

0.52 

89.1 

319 

84 

-1.82 

2.51 

0.217 

0.61 

133.0 

320 

84 

-0.60 

2.11 

0.142 

0.73 

221.0 

321 

84 

-0.62 

2.00 

0.123 

0.76 

263.0 

+  All  regressions  are  significant  at  P  <0.0001,  except  where  indicated: 
N/S  =  not  significant,  *  =  P  <0.05,  **  =  P  <0.01,  ***  -  P  <0.001. 


Table  2. - -Regression   analyses   of   the    log(variance)    against    log(mean)    of   tree 
ring  increments  before  and  after  19A3+ 


Site 
No. 


Pre/ 
post 


Log(a) 


Standard    Adjust   F- 


error  of  (b) 


ratio   Signif. 


Tennessee  Valley  Authority  Plots 


1 

0 

93 

1.04 

1.26 

0.158 

0.40 

63.2 

1 

42 

-0.43 

1.77 

0.546 

0.19 

10.6 

** 

3 

0 

46 

2.32 

0.60 

0.274 

0.08 

4.8 

* 

1 

42 

0.38 

1.47 

0.187 

0.60 

62.1 

4 

0 

50 

-0.65 

1.97 

0.341 

0.40 

33.4 

1 

42 

-2.28 

2.73 

0.190 

0.83 

206.0 

6 

0 

54 

3.48 

0.12 

0.300 

0 

0.2 

N/S 

1 

42 

0 

1.51 

0.535 

0.14 

8.0 

** 

7 

0 

66 

-0.65 

2.05 

0.220 

0.68 

86.5 

1 

42 

-0.65 

2.05 

0.220 

0.68 

86.5 

9 

0 

93 

-0.23 

1.79 

0.147 

0.61 

147.0 

1 

42 

-0.54 

2.00 

0.202 

0.70 

97.6 

10 

0 

93 

0.36 

1.44 

0.213 

0.33 

45.7 

1 

42 

4.35 

0 

0 

0 

1.0 

N/S 

12 

0 

83 

0.44 

1.49 

0.184 

0.44 

66.2 

1 

42 

-0.82 

2.04 

0.508 

0.27 

16.1 

*** 

17 

0 

38 

-0.91 

2.02 

0.560 

0.24 

13.0 

*** 

1 

42 

-1.03 

2.18 

0.226 

0.69 

93.0 

19 

0 

86 

0.08 

1.54 

0.160 

0.52 

92.0 

1 

42 

-0.90 

2.02 

0.225 

0.48 

38.7 

20 

0 

73 

0.14 

1.54 

0.171 

0.52 

79.6 

1 

42 

0.78 

1.26 

0.252 

0.37 

25.0 

22 

0 

58 

1.69 

0.87 

0.302 

0.11 

8.4 

** 

1 

42 

0.87 

1.28 

0.247 

0.39 

27.1 

23 

0 

93 

-0.96 

2.07 

0.178 

0.59 

134.0 

1 

42 

-0.91 

2.07 

0.292 

0.66 

80.4 

54 


Table  2 .- -Regression  analyses   of   the   log(variance)    against    log(mean)    of   tree 
ring  increments   before  and  after  1 943+- -Continued 


Site 

Pre/ 

Standard 

Adjust 

F- 

No. 

post 

N 

Log(a) 

b 

error  of  (b) 

R2 

ratio 

Signif . 

24 

0 

64 

-1.81 

2.46 

0.158 

0.79 

242.0 

1 

42 

-2.25 

2.66 

0.330 

0.61 

64.6 

25 

0 

93 

0.48 

1.42 

0.215 

0.32 

41.6 

1 

42 

0.79 

1.30 

0.204 

0.49 

40.7 

26 

0 

73 

1.76 

0.84 

0.188 

0.20 

19.5 

1 

42 

0.84 

1.23 

0.375 

0.19 

10.8 

** 

27 

0 

93 

1.05 

1.16 

0.350 

0.10 

11.1 

** 

1 

42 

-0.23 

1.83 

0.258 

0.55 

50.2 

28 

0 

90 

1.59 

0.84 

0.382 

0.04 

4.9 

* 

1 

42 

-2.68 

2.91 

0.443 

0.51 

43.1 

29 

0 

53 

-0.41 

1.85 

0.109 

0.64 

95.3 

1 

42 

-2.21 

3.04 

0.147 

0.91 

425.0 

34 

0 

83 

1.24 

1.10 

0.096 

0.62 

132.0 

1 

42 

-0.67 

1.77 

0.499 

0.22 

12.6 

**•*• 

36 

0 

93 

-3.31 

3.17 

0.370 

0.44 

73.4 

1 

42 

-4.67 

3.95 

0.181 

0.92 

477.0 

37 

0 

92 

-1.68 

0.79 

0.139 

0.25 

32.1 

Tennessee  Valley  Authority  Plots 


1 

42 

-1.19 

2.22 

0.253 

0.65 

76.5 

46 

0 

93 

-0.09 

1.79 

0.110 

0.74 

267.0 

1 

42 

-1.64 

2.48 

0.339 

0.56 

53.3 

54 

0 

93 

-0.66 

1.94 

0.471 

0.15 

16.9 

1 

42 

-2.66 

2.80 

0.511 

0.42 

30.1 

55 

0 

43 

-2.37 

2.71 

0.270 

0.70 

101.0 

1 

42 

-4.09 

3.57 

0.572 

0.48 

38.9 

132 

0 

48 

1.08 

1.25 

0.216 

0.36 

27.4 

1 

42 

-0.57 

1.85 

0.454 

0.28 

16 . 7   *** 

153 

0 

59 

2.28 

0.49 

0.231 

0.06 

4.4  * 

1 

42 

1.32 

1.14 

0.286 

0.27 

15 . 9  *** 

158 

0 

45 

1.48 

0.93 

0.108 

0.62 

74.0 

1 

42 

-2.95 

2.89 

0.370 

0.59 

60.9 

National  Park  Service  Plots 


301 

0 

42 

-1.50 

2.57 

0.421 

0.47 

37.2 

1 

42 

-1.96 

2.86 

0.152 

0.90 

353.0 

302 

0 

42 

-2.77 

2.86 

0.404 

0.54 

50.1 

1 

42 

-1.46 

2.40 

0.193 

0.79 

153.0 

303 

0 

42 

-1.40 

2.44 

0.236 

0.72 

108.0 

1 

42 

-1.50 

2.59 

0.211 

0.78 

150.0 

304 

0 

42 

-0.92 

2.11 

0.439 

0.35 

23.1 

1 

42 

2.58 

0.51 

0.261 

0.07 

3.9  * 

305 

0 

42 

-0.78 

2.04 

0.239 

0.64 

73.0 

1 

42 

0.29 

1.49 

0.159 

0.68 

87.2 

307 

0 

42 

-1.73 

2.61 

0.109 

0.93 

571.0 

1 

42 

0.27 

1.65 

0.194 

0.63 

72.3 

308 

0 

42 

-1.79 

2.58 

0.441 

0.45 

31.2 

1 

42 

-2.49 

3.06 

0.504 

0.47 

36.9 

309 

0 

42 

-3.51 

3.34 

0.558 

0.46 

35.8 

1 

42 

-1.19 

3.31 

0.272 

0.63 

72.1 

310 

0 

42 

-0.41 

1.92 

0.158 

0.78 

148.0 

1 

42 

-2.02 

2.77 

0.372 

0.57 

55.5 

311 

0 

42 

-2.43 

2.70 

0.224 

0.78 

145.0 

1 

42 

-0.18 

1.83 

0.166 

0.75 

122.0 

312 

0 

42 

2.40 

0.59 

0.261 

0.09 

5.1   * 

1 

42 

1.63 

0.84 

0.237 

0.22 

12.5   ** 

313 

0 

41 

-3.64 

3.16 

0.448 

0.53 

49.7 

1 

42 

5.19 

0.63 

0.585 

0 

1.2  N/: 

314 

0 

41 

0.30 

1.46 

0.333 

0.31 

19.1 

55 


Table  2 . - -Regression   analyses    of    the    log(variance)    against    log(mean)    of    tree 
ring  increments   before  and  after  1 943+- -Continued 


Site 
No. 


Pre/ 
post 


Log(a) 


Standard 
error  of  (b) 


Adjust  F- 


Rz 


ratio  Signif. 


National  Parks  Service  Plots 


315 


316 


317 


318 


319 


320 


321 


42 
42 
42 
42 
42 
42 
42 
42 
42 
42 
42 
42 
42 
42 
42 


2 

.35 

0 

.57 

4 

.42 

3 

.34 

1 

.14 

2 

.22 

4 

.28 

0 

.51 

0 

,95 

1 

.14 

0 

.18 

1 

,73 

0 

.39 

1. 

.87 

0 

.60 

1, 

.98 

0 

.09 

1. 

90 

2 

.14 

2. 

69 

2 

.99 

2, 

99 

1, 

.62 

2. 

61 

0 

86 

1. 

40 

0 

95 

2. 

17 

0 

12 

1, 

60 

0.357 
0.341 
0.207 
0.673 
0.249 
0.225 
0.281 
0.330 
0.209 
0.265 
0.380 
0.181 
0.183 
0.173 
0.191 


0 

.04 

2 

6 

0 

.72 

108 

0 

0 

.74 

115 

,0 

0 

0 

.6 

0 

.33 

21 

,1 

0 

.59 

59 

,5 

0 

.52 

44 

6 

0 

46 

35 

,8 

0 

67 

33 

,1 

0, 

,71 

103 

,0 

0 

60 

62 

.1 

0 

83 

208 

0 

0 

59 

59 

.1 

0 

79 

158 

0 

0 

63 

70 

1 

N/S 


N/S 


+A11  regressions  are  significant  at  P  <0.0001, except  where  indicated: 
N/S  =  not  significant,  *  =  P  <0.05,  **  =  F  <0.01,  ***  =  P  <0.001 
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Red  Spruce  Tree  Ring  Analysis  Using  a  Kalman  Filter 
Paul  C.  Van  Deusen 


SUMMARY 

A  Kalman  filter  was  applied  to  red  spruce 
(Picea  rubens  Sarg.)  tree  ring  data  collected 
from  the  Great  Smoky  Mountains  by  the  Tennessee 
Valley  Authority  and  the  National  Park  Service . 
A  new  standardization  method  was  developed  that 
can  be  justified  with  a  model-based  assumption. 
The  variance  of  the  standardized  growth 
chronology  appears  to  have  increased  in  recent 
years.  The  sensitivity  of  red  spruce  to  climate 
began  to  increase  in  the  late  1960's  and  leveled 
off  in  the  early  1980' s.  It  is  possible  that 
increasing  climatic  sensitivity  and  variance  are 
related  to  balsam  wooly  adelgid  activity  in  these 
stands . 

ISTRODUCTION 

Tree  ring  data  provide  one  of  the  few 
historical  records  for  scientists  to  assess  the 
impact  of  atmospheric  deposition  influences  on 
forests  (ADIF) .  However,  to  adequately  assess 
this  impact,  historical  information  on  weather 
and  pollution  levels  are  also  needed.  Although 
past  weather  records  are  available  from  weather 
stations,  deposition  levels  can  only  be  inferred 
from  proxy  variables  such  as  coal  consumption  of 
nearby  power  plants.  My  analysis  will  be  limited 
to  examining  the  trends  in  the  tree  ring  series 
and  attempting  to  explain  these  trends  with 
average  monthly  temperature  and  total 
precipitation  records.  The  Kalman  filter  is 
shown  to  be  a  useful  tool  for  this  type  of 
analysis. 

DATA 

The  National  Park  Service  (NPS)  and  the 
Tennessee  Valley  Authority  (TVA)  had  collected 
data  from  plots  located  in  the  Great  Smoky 
Mountains.  Spatial  relationships,  current 
diameter,  and  species  were  recorded  for  each  tree 
on  each  plot.  A  few  dominant  trees  were  selected 
at  each  NPS  plot  from  which  two  increment  cores 
were  taken,  but  pith  date  was  not  recorded;  tree 
rings  were  dated  back  to  1900.  Dominant  trees  on 
TVA  plots  were  not  cored,  but  pith  date  was 
recorded;  unfortunately  tree  ring  widths  were 
only  available  back  to  1850.  Information  from 
approximately  200  red  spruce  trees  were  available 
in  both  data  sets.  A  data  set  from  a  site  on 
Clingman's  Dome  (North  Carolina)  consisting  of  38 
cores  on  19  trees  was  also  utilized  in  this 
study.  Ed  Cook  of  the  Lamont-Doherty  Tree  Ring 
Laboratory  collected  and  cross -dated  this  data 
set.  The  TVA  and  NPS  tree  ring  data  were 
processed  at  Oak  Ridge  National  Laboratories. 

Generally,  only  data  that  would  be  commonly 
available  in  a  dendrochronological  study  were 
used;  these  data  included  ring  width,  elevation, 


pith  date,  and  regional  weather  data.  National 
Weather  Service  Climatic  Division  data  from 
stations  in  the  northern  mountains  of  North 
Carolina  were  summarized  to  provide  total 
rainfall  and  average  temperature  by  month 
beginning  in  1931. 

STANDARDIZATION 

A  basic  goal  of  many  tree  ring  studies  is  to 
produce  a  common  chronology  from  a  group  of  trees 
to  be  representative  of  a  site.  It  is  probable 
that  this  chronology  will  demonstrate  a  common 
signal  to  which  all  trees  in  the  area  have 
responded.  In  order  to  amplify  the  common 
signal,  an  attempt  is  made  to  eliminate 
individual  tree  signals  that  are  unrelated  to  the 
common  signal.  The  age-related  biological  growth 
signal  can  be  removed  by  a  number  of  procedures 
categorically  referred  to  as  standardization. 
Graybill  (1982)  presents  methods  where  a  growth 
model  is  individually  fit  to  each  tree  ring 
series. 

Graybill  (1982)  presents  a  hypothetical 
breakdown  of  the  raw  ring  width  for  a  single  tree 
at  time  t,  R(t) ,  as  follows: 


R(t) 


+  Bf-  +  Dl,-  +  D2r  +  e, 


(1) 


where 


is  unique 


Ct  is  the  macroclimatic  signal  common  to 
all  trees; 

Bt  is  the  biological  growth  trend  -  a 
function  of  tree  age ; 

Dlt  is  a  disturbance  signal  that 
to  the  individual ; 

D2t  is  a  disturbance  signal  common  to  most 
or  all  individuals- -possibly  caused  by 
fire,  insects,  or  pollution;  and 

et  accounts  for  random  disturbance. 


In  order  to  maximize  the  macrosignals  (C  and 
D2),  it  is  necessary  to  recognize  and  remove  the 
microsignals  (B  and  Dl)  as  much  as  possible  (Cook 
1987).   An  index  is  formed  as  followes: 


Index(t)  =  R(t)/Y(t), 


(2) 


where  Y(t)  is  the  model-based  prediction  of  R(t) . 
This  produces  a  new  index  series  with  an 
expectation  of  1,  a  more  homogeneous  variance, 
and  a  smaller  first  order  autocorrelation  than 
the  original  series  (Fritts  1976).  Graybill 
(1982)  presents  negative  exponentials  and 
orthogonal  polynomials  as  potential  prediction 
models. 

Cook  (1987)  discusses  the  use  of  splines  to 
replace  Graybill's  models.  There  is  a 
possibility  that  the  disturbance  signal  (Dl)  as 
well  as  the  B-  signal  may  be  removed  with  the 
spline  approach,  and  more  user  interaction  and 
expert  opinion  are   required.    Warren   (1980) 
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presents  an  alternate  model -based  approach  that 
could  also  be  used  to  remove  the  Dl  signal,  but 
it  also  requires  much  user  interaction. 


A  NEW  STANDARDIZATION  PROCEDURE 

A  method  was  sought  for  this  study  that 
required  little  subjectivity  and  could  be  used  to 
automatically  process  many  hundreds  of  cores . 
The  method  begins  with  a  standard  sigmoidal  model 
for  diameter  at  some  age,  D(A) : 


D(A)  =  b(l 


e-**) 


(3) 


where  b  is  the  asymptote  parameter,  and 

k  is  the  shape  parameter. 

Differentiating  model  (3)  with  respect  to  age 

gives  a  diameter  growth  model: 


d_PJA)_  .  b<rkA 
dA 

=  2  R(A) 


(4) 


Model  (4)  is  appropriate  for  radial  increment 
data  and  is  similar  to  standard 
dendrochronological  methods . 

Assuming  model  (4)  ,  two  steps  are  required  to 
remove  the  age-related  trend  from  the  data. 
First  take  the  natural  log  of  model  (4),  giving: 


log[R(A)]  =  constant  -kA  . 


(5) 


Then  take  the  first  differences  of  model  (5), 
giving : 


log[R(A)]  -  log  [R(A-l)]  =  -k 

=  Alog[R(t)] 


(6) 


Thus  a  simple  transformation  removes  the  age- 
related  trend  from  the  tree  ring  series .  Result 
(6)  can  be  justified  intuitively  by  viewing  it  as 
a  relative  growth  rate.  Taking  the  log  of  R(t) 
puts  it  on  a  relative  scale,  and  the  first 
difference  is  just  a  numerical  first  derivative 
that  can  be  viewed  as  a  growth  rate .  This 
transformation  can  be  quickly  performed  without 
sophisticated  software  or  user  interaction,  an 
important  advantage  for  large  data  sets . 

Plotting  the  data  before  and  after 
transformation  indicated  that  the  new  series  was 
stationary  and  that  the  age-related  trend  had 
been  removed  as  expected  (fig.  1).  Figures  1  and 
2  show  how  the  transformation  creates  similar 
series  from  a  young  tree  and  an  old  tree  that 
looked  quite  different  before  transformation. 
Notice  that  1937,  1969,  and  1981  are  all  low  on 
the  transformed  series. 

The  analysis  can  then  proceed  on  the 
transformed  data  by  assuming  that  the  new  series 
is  composed  of  the  following  signals: 


Alog[R(t)]  =  Ct  +  Dlt  +  D2t  +  et, 


(7> 


where  Ct,  Dlt,  D2t,  and  et  are  defined  as  ir 
model  (1) . 

The  age-related  signal  is  now  removed.  Two 
macros ignal  terms  that  are  of  interest  remain  (C 
and  D2) ,  as  well  as  two  uninteresting  micros ignal 
terms  (Dl  and  e)  that  will  be  treated  as  noise. 


KALMAN  FILTER 

A  system  for  updating  and  predicting  is 
presented  by  Kalman  (1960)  based  on  a  linear 
dynamic  model.  These  models  are  generalizations 
that  can  generate  any  of  the  class  of  ARMA  models 
(Box  and  Jenkins  1976) ,  standard  multiple 
regression  models,  and  regression  models  with 
time  varying  parameters  (Harvey  1981) .  Applied 
to  dendrochronology,  the  Kalman  filter  provides  a 
means  of  simultaneously  reducing  a  number  of 
series  to  a  single  chronology  and  generating 
climate -based  predictions.  Furthermore,  the 
climate  parameters  can  be  allowed  to  vary  over 
time  to  provide  a  test  of  the  uniformitarian 
principle  that  conditions  in  the  past  are  similar 
to  the  present.  The  uniformitarian  principle  is 
the  fundamental  justification  for  the  use  of 
dendrochronology  to  infer  past  conditions . 

The  Kalman  filter  can  be  derived  from  Bayesian 
theory  (Harrison  and  Stevens  1976;  Meinhold  and 
Singpurwalla  1983)  or  with  least  squares  methods 
presented  by  Duncan  and  Horn  (1972) .  The 
equations  needed  to  implement  the  Kalman  filter 
are  presented  below,  and  the  reader  is  referred 
to  the  citations  for  the  theoretical  development. 

The  basic  difference  between  the  Kalman  filter 
and  usual  regression  models  is  that  the 
parameters  are  allowed  to  vary  over  time .  The 
relationship  between  the  vector  of  observed 
standardized  ring  widths  at  time  t  (Yt)  and  the 
parameters  (at)  *-s  called  the  observation  or 
measurement  equation: 

Yt  =  Ftat  +  vt, 

where  the  matrix  Ft  is  fixed  and  of  order  ntxp , 
at   is   a   pxl   vector   of   underlying   state 
parameters,  and  vt  is  an  ntxl  vector  of  residuals 
with  zero  expectation  and  variance  matrix  Vt. 


The   state  variable , 


evolves  over  time 


according  to  a  first  order  Markov  process  as 
defined  by  the  transition  or  system  equation: 

at  =  Gtat-1  +  wt  , 

where  Gt  is  a  fixed  pxp  matrix,  and  wt  is  a  pxl 
vector  of  residuals  with  zero  expectation  and 
variance  matrix  Wt. 

The  error  terms  vt  and  wt  are  assumed  to  be 
independent  white  noise  series.  The  quantities 
that  must  be  known  include  the  matrices  (Ft  and 
Gt)  that  premultiply  the  state  variables  (at  and 
at-l)  .  The  matrices  Ft  and  Gt  correspond  to 
independent  variables  in  regression  theory,  or, 
when  dealing  with  rocket  trajectories,  they  come 
from  well-defined  physical  laws.  The  more 
complex  problem  comes  from  the  need  to  know  the 
variance  matrices  (Wt  and  Vt)  ,  because  in  many 
statistical  applications  this  will  require  some 
user  subjectivity. 

The  equations  needed  to  estimate  the  state 
variables  can  be  divided  into  three  parts : 
prediction  equations,  updating  equations,  and 
smoothing  equations.  Let  at.^  denote  the  optimal 
estimator  of  at-i  based  on  all  information  up  to 
and  including  Yt-i-  The  covariance  matrix  of 
at-i  -  <*t-l  will  be  Pt-1-  The  prediction 
equations  for  at  and  the  associated  covariance 
matrix  given  at-i  and  Pt-1  are: 
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Figure  I. --Ring  widths  of  a  relatively  young  red  spruce  in  raw  (upper  graph) 
and  standardized  (lower  graph)  form  for  a  young  tree.  First 
differences  of  the  natural  log  were  used  to  produce  the  lower 
standardized  series. 


at/t-l  "  Gtat-1   -  and 
pt/t-l  "  GtPt.iGt'  +  Wt 


(3a) 
(3b) 


When  Yt  becomes  available,  the  updating 
equations  for  the  estimate  of  at  and  the 
associated  covariance  matrix  are: 


-1 


Et ,  and 
-  pt/t-l  -  Pt/t-lFt'Zt-J-Ft-P,-/,-.! 
where 


"  at/t-l  +  Pt/t-lFt'Zt, 


(4a) 
(4b) 


Et  -  Yt  -  Ftat/t.x, 
Zt  =  FtPt/t.^Ft'  +  Vt  ,  and,  for  computational 
savings , 


-1 


V 


■r1    -   Vr;1Ff[F^Vf--1Ft 


't,  rtlrt  vt,  rt 

■l]-IFt»l 

The  estimate  of  at  in  (4a)  is  the  sum  of  its 


+  Pt/t-l"J-]"xFt'Vt-1. 


estimate  at  time  t-1  and  a  weighted  average  of 
the  prediction  errors  (Et)  . 

At  any  time,  at  is  an  optimal  estimate,  given 
all  previous  information,  but  only  the  estimate 
at  time  T  (the  final  period)  contains  all  of  the 
information  available.  Given  all  information, 
the  optimal  solution  for  any  time  t  is  referred 
to  as  signal  extraction  or  smoothing.  Smoothing 
begins  with  the  solutions  at  time  T  and 
recursively  goes  backwards  to  time  1.  This 
yields  the  optimal  smoothed  estimates  of  the 
state  variables  with  associated  covariance 
matrices  as  follows: 


»t/T  ~  at  +  pt  (at+l/T  "  Gt+lat).  and 


t/T 


-  P 


t  +  *t  <Pt+l/T  "  pt+l/t)Pt 


*, 


(5a) 
(5b) 
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Figure  2 . - -Ring  widths  of  an  older  red  spruce  in  raw  (upper  graph)  and 
standardized  (lower  graph)  form  for  an  old  tree.  First 
differences  of  the  natural  log  were  used  to  produce  the  lower 
standardized  series. 


where 

aT/T  =  aT  f°r  t*le  starting  value  on  the 

state  variables, 
PT/T  =  PT  f°r  t*le  covariance  starting 

values  and 

Pt*  "  PtGt+l'  Pt+l/t"1 
APPLYING  THE  KALMAN  FILTER  IN  DENDROCHRONOLOGY 

The  Kalman  filter  provides  a  complete  system 
for  prediction,  updating,  and  smoothing  that 
should  appeal  to  the  dendrochronologist .  In 
particular,  the  usual  procedure  of  climate 
prediction  from  a  single  chronology  formed  from 
many  standardized  tree  ring  series  could  be 
refined.  Two  applications  of  the  Kalman  filter 
will  be  presented.     The   first  shows  how  a 


chronology  can  be  formed  simultaneously  with  its 
climate  predictor,  and  the  second  is  an 
application  to  dendroclimatology.  For  an 
additional  application,  readers  should  see  Jones 
(1980)  where  various  models  are  fit  to  drought 
data  reconstructed  from  tree  rings . 

Application  1 

The  Kalman  filter  was  applied  to  the  red  spruce 
data  collected  by  the  NPS  and  TVA.  The  objective 
was  to  simultaneously  form  a  single  chronology 
and  its  climate -based  prediction.  The  data  were 
first  standardized  by  taking  the  first  difference 
of  the  logarithms .  The  following  Kalman  filter 
was  then  applied: 


60 


Observation  equation  Yt  -  F, 


[:::] 


+  v. 


Transition  equation  a^t  =  CtQ2  t-1  +  wlt 
a2t  =   a2,t-l  +  w2t 


(6a) 

(6b) 
(6c) 


of  observed  tree  ring 


where 

Yt  is  the 

data  at  time  t, 

Ft  is  a  matrix  with  two  columns  of  length 

nt  with  the  first  column  being  l's  and 

the  second  O's, 

vt>  wlt>  an<*  w2t  are  random  errors, 

Ct   is  the  climate  variable, 

a^t  is  the  value  of  the  chronology  at  time 

t ,  and 

a2t  is  the  climate  effect  at  time  t. 


The 


aKi 


matrix  Vt  was  defined  to 
estimated  from  vector  Yt 


be 
as 


covariance 
£j.£  where  °    t.    was  eaLj.uia.i_eu  liuu  veti-ut  i£ 

S(Y£t  -Yt)/[nt-l]  ,  and  It  was  an  identity  matrix 

of  order  nt.    In  other  words,  the  trees  were 

treated  a  priori  as  independent.  The  covariance 
matrix  Wt  was  defined  as 


*-*[V..kO 


Thus  the  state  variables  were  assumed  to  be  a 
priori  independent.  The  elements  on  the  diagonal 
were  chosen  to  allow  the  state  variables  to  vary 
enough  to  respond  to  a  trend,  but  not  enough  to 
absorb  random  fluctuations.  Since  a^t  might 
approximate  the  average  of  the  vector  Yt,  n^ 
was  chosen  for  the  upper  diagonal,  and  the  factor 
0.01  was  used  in  the  lower  diagonal  to  prevent 
c*2t  from  fluctuating  in  response  to  random 
disturbances.  The  results  were  robust  to  changes 
in  the  W  and  V  matrices,  which  lends  credence  to 
these  values . 

Data  from  NPS  plots  dates  back  to  1900.  The 
above  algorithm  was  applied  to  the  data  from  1901 
through  1983,  because  taking  first  differences 
eliminates  the  first  observation.  The  climate 
data  were  available  from  1931  through  1983,  so 
the  filter  was  modified  during  the  "preclimate" 
period  (1901-30)  by  setting  W2t  and  Ct  to  zero. 

Results  of  Application  1 

Climate  lagged  by  1  year  was  found  to  best 
predict  the  chronology  or  time  trend  in  the  data. 
The  climate  variable  used  was  a  linear 
combination  of  all  the  monthly  rainfall  and 
temperature  data  for  a  year,  as  described  in 
table  1 . 

Figure  3  presents  the  chronologies  for  the 
three  data  sets.  The  NPS  and  Clingman's  Dome 
chronologies  were  very  similar,  particularly  over 
the  past  10  years.  Variation  in  the  standardized 
series  had  a  tendency  to  increase  in  the  more 
recent  years.  This  increase  was  most  evident  in 
the  TVA  chronology,  but  there  was  no  suggestion 
for  the  cause  of  increases. 

The  time  trend  predictions  based  on  lagged 
climate  and  the  climate  effect  variable  c*2t  with 
95  percent  confidence  intervals  are  plotted  in 
figures  4,  5,  and  6.  Lagged  climate  predicted 
the   chronologies   accurately  beginning   in   the 


Table  1.  The  climate  variable  used  was  a  linear 
combination  of  monthly  rainfall  and  temperature 
averages.  The      weight      below      is      the      value 

multiplied  by  the  corresponding  monthly  average 
to  create  this  combined  climate  variable .  The 
principle  components  method  was  used  to  create 
this  linear  combination.  September  temperature 
and  October  rainfall  have   the  largest  weights. 


Variable 


Month 


Weight 


Temperature 


Rainfall 


J 
F 
M 
A 
M 
J 
J 
A 
S 
0 
N 
D 

J 
F 
M 
A 
M 
J 
J 
A 
S 
0 
N 
D 


0.056 

-0.176 

-0.124 

-0.257 

-0.106 

-0.112 

-0.270 

-0.264 

-0.409 

0.013 

0.204 

0.040 

-0.056 
0.164 
0.083 

-0.039 
0.161 
0.199 
0.061 
0.026 
0.027 

-0.517 
0.166 
0.179 


1960's,  but  did  poorly  before  this.  The 
corresponding  plot  of  the  parameter  (the  climate 
effect)  that  multiplies  climate  also  showed  an 
increase  over  time.  This  indicated  that  the 
trees  were  becoming  more  responsive  to  climate  in 
the  1960 's  than  they  were  previously.  The  reason 
for  the  increased  sensitivity  to  climate  cannot 
be  determined  from  this  study.  One  might 
speculate  that  this  was  caused  by  thinning  in  the 
stands  as  a  result  of  insect  damage  to  the  fir 
component.  whether  this  is  related  to  pollution 
has  not  been  determined. 

Application  2 

The  Kalman  filter  can  be  applied  to 
simultaneous  prediction  of  past  climate  and  the 
single  common  chronology.  This  traditionally 
involves  averaging  the  individual  series  together 
as  a  first  step  to  form  the  single  chronology. 
The  climate  prediction  is  then  a  separate  step 
that  takes  place  without  the  complete  information 
contained  in  the  original  series. 

A  Kalman  filter  can  be  formulated  to  handle 
these  steps  simultaneously.  This  incorporates 
the  full  information  contained  in  the  data  while 
automatically  providing  a  prediction  system  for 
past  climate  with  associated  prediction 
intervals.    Although  the  simple  solution  given 
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Figure  3 . --Standardized  chronologies  from  the  Clingman's  Dome,   National  Park 
Service,      and     Tennessee     Valley     Authority     data     sets.  First 

differences  of  the  log  transform  were  used  to  standardize  each 
tree  ring  series.  A  Kalman  filter  was  used  to  produce  the 
chronology  as  described  in  equations    (6a) ,    (6b) ,    and   (6c) . 
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Figure  4. --Time  trend  predictions:  A,  Clingman's  dome  chronology  (solid 
line)  and  its  climate  based  prediction  (dashed  line)  using  the 
Kalman  filter  described  by  equations  (6a),  (6b),  and  (6c);  B,  the 
trend  in  the  climate  parameter  given  in  equation  (6a) .  The 
climate  variable  is  a  principle  component  (linear  combination)  of 
monthly  average  rainfall   and   temperature  variables . 
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Figure  5. --Time  trend  predictions:  A,  National  Park  Service  chronology 
(solid  line)  and  its  climate  based  prediction  (dashed  line)  using 
the  Kalman  filter  described  by  equations  (6a),  (6b),  and  (6c);  B, 
the  trend  in  the  climate  parameter  with  95  percent  confidence 
intervals,  given  in  equation  (6a).  The  climate  variable  is  a 
principle  component  (linear  combination)  of  monthly  average 
rainfall  and  temperature  variables. 
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Figure  6. --Time  trend  predictions:  A,  Tennessee  Valley  Authority  chronology 
(solid  line)  and  its  climate  based  prediction  (dashed  line)  using 
the  Kalman  filter  described  by  equations  (6a),  (6b),  and  (6c);  B, 
the  trend  in  the  climate  parameter  given  in  equation  6a.  The 
climate  variable  is  a  principle  component  (linear  combination)  of 
monthly  average   rainfall    and    temperature   variables . 
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below  does  not  achieve  the  potential  of  the 
method,  it  was  chosen  to  emphasize  some  important 
points.   Consider  the  following  formulation: 


Observation 
equations 


Transition 
equations 


fit"!  ro  »ir*»-'i  r""i  "c> 
L«,J-Lo  iJL«2.wJl-»J  (7d) 


where 


Yt  is  the  ntxl  vector  of  standardized  ring 

widths  at  time  t, 
Ct  is  a  qxl  vector  of  observed  climate 

variables , 
jt  is  a  vector  of  l's  of  length  nt , 
0t  is  a  vector  of  O's  of  length  nt , 


<*!(-  is  the  value  of  the  chronology  at  time 

t, 
c*2t  is  tne  climate  parameter  at  time  t,  and 
vlt>  v2t»  wlt>  w2t  are  random  errors. 

Equations  (7a)  through  (7d)  define  a  system  that 
handles  the  dendroclimatological  problems  of 
forming  a  single  chronology  for  a  site  providing 
a  means  of  predicting  past  climate,  and  testing 
the  uniformitarian  principle. 

Some  details  on  implementing  the  above  system 
must  be  noted.  The  iterative  process  is  started 
at  time  T  rather  than  time  1,  since  past  climate 
predictions  are  needed.  Suppose  that  the  climate 
variables  are  available  from  time  T  to  time  t  , 
where  1  <  t*  <  T.  During  this  interval  of  known 
climate  (INTC)  the  parameter  c*2  is  estimated. 
Beyond  INTC ,  equation  (7b)  is  eliminated  and  the 
Kalman  filter  is  allowed  to  generate  new  values 
of  c*2  back  to  time  1  (the  earliest  time  for  which 
tree  data  are  available)  and  simultaneously 
produce  the  chronology  and  climate  predictions  as 
a2t^t+l-  Furthermore,  the  trend  can  be  examined 
in  the  c*2  parameter  over  INTC  to  see  if  the 
uniformitarian  principle  is  valid,  although  the 
Kalman  filter  will  tend  to  move  02  along  the 
established  trend  making  the  uniformitarian 
assumption  less  important. 


Results  of  Application  2 

Starting  values  must  be  supplied  to  the  system  in 
(7a)  through  (7d)  for  the  parameter  vector  and 
the  variance  matrices  Vt  and  Wt;  Vt  was  assumed 
diagonal  with  variances  estimated  from  the 
vectors  Yt  as  in  example  1.  The  last  diagonal 
element  in  Vt  is  the  variance  of  V2t  and  was 
estimated  from  the  variance  in  the  known  climate 
data.  The  matrix  Wt  was  also  assumed  diagonal 
and  chosen  similarly  to  example  1. 

A  principle  component  of  climate  variables  was 
used  for  Ct.  A  starting  value  of  zero  was  used 
for  the  chronology  parameter  a^ ,  and  «2  was 
started   at   10.     Unfortunately,   climate   was 


predicted  poorly  (fig.  7).  The  values  before  1931 
were  predictions  generated  by  the  model.  Figure 
7  shows  the  climate  parameter  c*2  trending  over 
time  with  a  95  percent  confidence  interval.  The 
parameter  was  tending  toward  zero  as  the 
confidence  interval  expanded,  which  is  not 
surprising  given  the  poor  predictions.  Although 
the  model  presented  in  equations  (7a)  through 
(7d)  may  not  be  ideal,  this  demonstrates  how  one 
could  use  the  Kalman  filter  for  predicting  past 
climate . 

CONCLUSION 

The  focus  of  this  paper  was  to  apply  the  Kalman 
filter  to  the  study  of  tree  rings.  The  Kalman 
filter  provides  a  natural  way  of  handling  many  of 
the  problems  that  dendrochronologists  encounter. 


CLINGMAN'S  DOME  DATA-PREDICTED  VS  ACTUAL  CLIM 
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Figure  7 . - -Climate  prediction:  A,  prediction 
(dashed  line)  of  past  climate  obtained 
from  tree  ring  data  using  the  Kalman 
filter  described  by  equations  (7a) 
through  (7d) .  The  solid  line  is  known 
climate,  which  is  available  back  to 
1931;  B,  the  climate  parameter  trend 
and  95  percent  confidence  interval. 
This  plot  indicates  that  the 
uniformitarian  principle  does  not  hold 
here,  since  the  parameter  is  tending 
toward  zero. 
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The  first  application  involved  prediction  of  an 
average  chronology  with  climate  incorporated  in 
the  process.  The  parameter  associated  with  the 
climate  was  allowed  to  vary  over  time.  The 
climate  parameter  followed  a  sigmoid  curve  that 
implied  an  increasing  sensitivity  to  climate  with 
time.  One  might  speculate  that  insect-caused 
thinning  of  the  fir  component  accounts  for  this 
phenomenon. 

The  second  application  gave  some  insight  as  to 
how  the  method  can  be  applied  to  more  typical 
dendrochronological  needs,  i.e.,  predicting  the 
past  value  of  climate  from  tree  rings.  It  was 
previously  indicated  in  application  1  that  the 
data  set  employed  was  insensitive  to  climate  in 
the  past  and  thus  predictions  were  poor. 
Although  the  filter  was  not  extremely 
sophisticated,  it  showed  how  one  might  proceed 
with  such  a  problem.  Future  mean  tree  ring 
values  were  used  to  predict  current  climate  while 
simultaneously  developing  the  mean  chronology. 
The  trend  in  the  climate  parameter  could  also  be 
inspected  as  an  indication  of  the  validity  of  the 
uniformitarian  principle. 
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Four  different  analyses  of  red  spruce  tree  ring  data 
from  the  Great  Smoky  Mountains  are  presented  along  with 
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